| 8 | #include <iostream> |
| 9 | |
| 10 | bool ODEIntegrators::AdaptiveRK54(AbstractODEIntegrator& ode, double tmin, double tmax, double hmin, double hmax, double eps_allowed, |
| 11 | double step_relax) { |
| 12 | // Get the starting array of variables of integration |
| 13 | std::vector<double> xold = ode.get_initial_array(); |
| 14 | const long N = static_cast<long>(xold.size()); |
| 15 | |
| 16 | // Start at an index of 0 |
| 17 | int Itheta = 0; |
| 18 | double t0 = tmin; |
| 19 | double h = hmin; |
| 20 | |
| 21 | // Figure out if t is increasing or decreasing in the integration and set a flag |
| 22 | bool forwards_integration = ((tmax - tmin) > 0); |
| 23 | // If backwards integration, flip the sign of the step |
| 24 | if (!forwards_integration) { |
| 25 | h *= -1; |
| 26 | } |
| 27 | |
| 28 | double max_error = NAN; |
| 29 | |
| 30 | 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); |
| 31 | |
| 32 | // t is the independent variable here, where t takes on values in the bounded range [tmin,tmax] |
| 33 | do { |
| 34 | |
| 35 | // Check for termination |
| 36 | bool abort = ode.premature_termination(); |
| 37 | if (abort) { |
| 38 | return abort; |
| 39 | } |
| 40 | |
| 41 | bool stepAccepted = false, disableAdaptive = false; |
| 42 | |
| 43 | while (!stepAccepted) { |
| 44 | |
| 45 | // reset the flag |
| 46 | disableAdaptive = false; |
| 47 | |
| 48 | // If the step would go beyond the end of the region of integration, |
| 49 | // just take a step to the end of the region of integration |
| 50 | if (forwards_integration && (t0 + h > tmax)) { |
| 51 | disableAdaptive = true; |
| 52 | h = tmax - t0; |
| 53 | } |
| 54 | if (!forwards_integration && (t0 + h < tmax)) { |
| 55 | disableAdaptive = true; |
| 56 | h = tmax - t0; |
| 57 | } |
| 58 | |
| 59 | ode.pre_step_callback(); |
| 60 | |
| 61 | // We check stepAccepted again because if the derived class |
| 62 | // sets the variable stepAccepted, we should not actually do the evaluation. |
| 63 | // cppcheck-suppress identicalInnerCondition |
| 64 | if (!stepAccepted) { |
| 65 | |
| 66 | Eigen::Map<Eigen::VectorXd> xold_w(&(xold[0]), N); |
| 67 |
nothing calls this directly
no test coverage detected