Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/tests/odeint/embeddedrungekutta.cpp
7178ca623984reactmicp/unsaturated

Last Change
Raw File

embeddedrungekutta.cpp
View Options

#include "catch.hpp"
#include "odeint/runge_kutta_step.hpp"
using namespace specmicp::odeint;
void rhs1(double x, const double& yin, double& yout)
{
yout = -2*yin;
}
void rhs2(double x, const vector_t<2>& y, vector_t<2>& yout)
{
yout(0) = y(1);
yout(1) = -y(0);
}
TEST_CASE("Embedded Runge Kutta")
{
SECTION("Simple test Cash-Karp")
{
CashKarpStep<1> algo(rhs1, butcher_cash_karp45);
double yout, yerr;
double tf = 0.01;
algo.run_step(1, -2, 0, tf, yout, yerr);
CHECK(yout == Approx(std::exp(-2*tf)).epsilon(1e-3));
tf = 0.05;
algo.run_step(1, -2, 0, tf, yout, yerr);
CHECK(yout == Approx(std::exp(-2*tf)).epsilon(1e-2));
}
SECTION("Simple test Dormand-Prince")
{
DormandPrinceStep<1> algo(rhs1, butcher_dormand_prince45);
double yout, yerr;
double tf = 0.005;
algo.run_step(1, -2, 0, tf, yout, yerr);
CHECK(yout == Approx(std::exp(-2*tf)).epsilon(1e-3));
tf = 0.01;
algo.run_step(1, -2, 0, tf, yout, yerr);
CHECK(yout == Approx(std::exp(-2*tf)).epsilon(1e-2));
}
SECTION("Spring Cash-Karp")
{
CashKarpStep<2> algo = get_cash_karp_step<2>(rhs2);
vector_t<2> yout, yerr;
vector_t<2> y0, dydx;
y0 << 1, 0;
dydx << 0, -1;
double tf=0.01;
algo.run_step(y0, dydx, 0, tf, yout, yerr);
CHECK(yout(0) == Approx(std::cos(tf)).epsilon(1e-3));
}
SECTION("Spring Dormand-Prince")
{
DormandPrinceStep<2> algo = get_dormand_prince_step<2>(rhs2);
vector_t<2> yout, yerr;
vector_t<2> y0, dydx;
y0 << 1, 0;
dydx << 0, -1;
double tf=0.01;
algo.run_step(y0, dydx, 0, tf, yout, yerr);
CHECK(yout(0) == Approx(std::cos(tf)).epsilon(1e-3));
}
SECTION("Adaptative test Cash-Karp")
{
CashKarpStep<1> algo = get_cash_karp_step<1>(rhs1);
vector_t<1> yout, dydx;
yout = 1;
dydx = -2;
StepLength stepl(0.5);
double x = 0.0;
int cnt = 0;
while (x < 5)
{
algo.rk_step(yout, dydx, x, stepl);
cnt +=1;
CHECK(yout == Approx(std::exp(-2*x)).epsilon(1e-3));
stepl.get_next();
algo.get_rhs(x, yout, dydx);
}
}
SECTION("Adaptative test Dormand-Prince")
{
DormandPrinceStep<1> algo = get_dormand_prince_step<1>(rhs1);
vector_t<1> yout, dydx;
yout = 1;
dydx = -2;
StepLength stepl(0.5);
double x = 0.0;
int cnt = 0;
while (x < 5)
{
algo.rk_step(yout, dydx, x, stepl);
cnt +=1;
CHECK(yout == Approx(std::exp(-2*x)).epsilon(1e-3));
stepl.get_next();
algo.get_rhs(x, yout, dydx);
}
}
SECTION("Cosine adaptative test - Cash-Karp")
{
CashKarpStep<2> algo = get_cash_karp_step<2>(rhs2);
vector_t<2> yout, dydx;
yout << 1, 0;
dydx << 0, -1;
StepLength stepl(0.5);
double x = 0.0;
int cnt = 0;
while (x <4)
{
algo.rk_step(yout, dydx, x, stepl);
cnt +=1;
CHECK(yout(0) == Approx(std::cos(x)).epsilon(1e-3));
stepl.get_next();
algo.get_rhs(x, yout, dydx);
}
}
}
c4science · Help