12 const long N =
static_cast<long>(xold.size());
20 bool forwards_integration = ((tend - tstart) > 0);
22 if (!forwards_integration) {
28 std::vector<double> xnew1(N), xnew2(N), xnew3(N), xnew4(N), xnew5(N), f1(N), f2(N), f3(N), f4(N), f5(N), f6(N), error(N), xnew(N);
39 bool stepAccepted =
false, disableAdaptive =
false;
41 while (!stepAccepted) {
44 disableAdaptive =
false;
48 if (forwards_integration && (t0 + h > tend)) {
49 disableAdaptive =
true;
52 if (!forwards_integration && (t0 + h < tend)) {
53 disableAdaptive =
true;
63 Eigen::Map<Eigen::VectorXd> xold_w(&(xold[0]), N);
65 if (std::abs(h) < hmin && !disableAdaptive) {
67 h = (forwards_integration) ? hmin : -hmin;
68 disableAdaptive =
true;
77 Eigen::Map<Eigen::VectorXd> xnew1_w(&(xnew1[0]), N), f1_w(&(f1[0]), N);
78 xnew1_w = xold_w + h * (1.0 / 5.0) * f1_w;
80 ode.
derivs(t0 + 1.0 / 5.0 * h, xnew1, f2);
81 Eigen::Map<Eigen::VectorXd> xnew2_w(&(xnew2[0]), N), f2_w(&(f2[0]), N);
82 xnew2_w = xold_w + h * (+3.0 / 40.0 * f1_w + 9.0 / 40.0 * f2_w);
84 ode.
derivs(t0 + 3.0 / 10.0 * h, xnew2, f3);
85 Eigen::Map<Eigen::VectorXd> xnew3_w(&(xnew3[0]), N), f3_w(&(f3[0]), N);
86 xnew3_w = xold_w + h * (3.0 / 10.0 * f1_w - 9.0 / 10.0 * f2_w + 6.0 / 5.0 * f3_w);
88 ode.
derivs(t0 + 3.0 / 5.0 * h, xnew3, f4);
89 Eigen::Map<Eigen::VectorXd> xnew4_w(&(xnew4[0]), N), f4_w(&(f4[0]), N);
90 xnew4_w = xold_w + h * (-11.0 / 54.0 * f1_w + 5.0 / 2.0 * f2_w - 70.0 / 27.0 * f3_w + 35.0 / 27.0 * f4_w);
92 ode.
derivs(t0 + h, xnew4, f5);
93 Eigen::Map<Eigen::VectorXd> xnew5_w(&(xnew5[0]), N), f5_w(&(f5[0]), N);
96 + h * (1631.0 / 55296 * f1_w + 175.0 / 512.0 * f2_w + 575.0 / 13824.0 * f3_w + 44275.0 / 110592.0 * f4_w + 253.0 / 4096.0 * f5_w);
99 ode.
derivs(t0 + 7.0 / 8.0 * h, xnew5, f6);
100 Eigen::Map<Eigen::VectorXd> xnew_w(&(xnew[0]), N), f6_w(&(f6[0]), N);
101 xnew_w = xold_w + h * (37.0 / 378.0 * f1_w + 250.0 / 621.0 * f3_w + 125.0 / 594.0 * f4_w + 512.0 / 1771.0 * f6_w);
103 Eigen::Map<Eigen::VectorXd> error_w(&(error[0]), N);
106 * (-277.0 / 64512.0 * f1_w + 6925.0 / 370944.0 * f3_w - 6925.0 / 202752.0 * f4_w - 277.0 / 14336.0 * f5_w + 277.0 / 7084.0 * f6_w);
108 max_error = error_w.norm();
112 if (disableAdaptive) {
117 if (max_error > eps_allowed) {
121 h *= std::min(step_relax * pow(eps_allowed / max_error, 0.3), 0.999);
122 stepAccepted =
false;
129 std::cout <<
format(
"accepted");
141 if (max_error < eps_allowed && disableAdaptive == false && max_error > 0) {
144 h *= step_relax * pow(eps_allowed / max_error, 0.2);
148 if (forwards_integration) {
149 h = std::min(h, hmax);
151 h = -std::min(std::abs(h), hmax);
155 if (forwards_integration && (t0 - tend > +1e-3)) {
158 if (!forwards_integration && (t0 - tend < -1e-3)) {
161 }
while (((forwards_integration) && t0 < tend - 1e-10) || ((!forwards_integration) && t0 > tend + 1e-10));
167 #if defined(ENABLE_CATCH)
168 # include <catch2/catch_all.hpp>
170 TEST_CASE(
"Integrate y'=y",
"[ODEIntegrator]") {
174 std::vector<double> t, h, y;
177 return std::vector<double>(1, 1);
185 this->t.push_back(t);
186 this->h.push_back(h);
187 this->y.push_back(y[0]);
194 virtual void derivs(
double t, std::vector<double>& y, std::vector<double>& yprime) {
199 SimpleODEIntegrator simple;
201 double yfinal_integration = simple.y[simple.y.size() - 1];
202 double tfinal_integration = simple.t[simple.t.size() - 1];
204 double yfinal_analytic = exp(4.0);
205 double error = yfinal_integration / yfinal_analytic - 1;
207 CAPTURE(yfinal_analytic);
208 CAPTURE(yfinal_integration);
209 CAPTURE(tfinal_integration);
210 CHECK(std::abs(error) < 1e-6);
211 CHECK(std::abs(tfinal_integration - 4) < 1e-10);