31 #ifndef _adevs_rk_45_h_
32 #define _adevs_rk_45_h_
33 #include "adevs_hybrid.h"
43 template <
typename X>
class rk_45:
55 double integrate(
double* q,
double h_lim);
56 void advance(
double* q,
double h);
66 double trial_step(
double h);
71 ode_solver<X>(sys),err_tol(err_tol),h_max(h_max),h_cur(h_max)
73 for (
int i = 0; i < 6; i++)
74 k[i] =
new double[sys->
numVars()];
75 dq =
new double[sys->
numVars()];
76 qq =
new double[sys->
numVars()];
85 for (
int i = 0; i < 6; i++)
93 while ((dt = integrate(q,h)) < h) h -= dt;
100 double err = DBL_MAX, h = std::min(h_cur*1.1,std::min(h_max,h_lim));
103 for (
int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
107 if (err <= err_tol) {
108 if (h_cur <= h_lim) h_cur = h;
113 double h_guess = 0.8*pow(err_tol*pow(h,4.0)/fabs(err),0.25);
114 if (h < h_guess) h *= 0.8;
119 for (
int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
123 template <
typename X>
127 this->sys->der_func(qq,dq);
128 for (
int j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
130 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.5*k[0][j];
131 this->sys->der_func(t,dq);
132 for (
int j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
134 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.25*(k[0][j]+k[1][j]);
135 this->sys->der_func(t,dq);
136 for (
int j = 0; j < this->sys->numVars(); j++) k[2][j] = step*dq[j];
138 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] - k[1][j] + 2.0*k[2][j];
139 this->sys->der_func(t,dq);
140 for (
int j = 0; j < this->sys->numVars(); j++) k[3][j] = step*dq[j];
142 for (
int j = 0; j < this->sys->numVars(); j++)
143 t[j] = qq[j] + (7.0/27.0)*k[0][j] + (10.0/27.0)*k[1][j] + (1.0/27.0)*k[3][j];
144 this->sys->der_func(t,dq);
145 for (
int j = 0; j < this->sys->numVars(); j++) k[4][j] = step*dq[j];
147 for (
int j = 0; j < this->sys->numVars(); j++)
148 t[j] = qq[j] + (28.0/625.0)*k[0][j] - 0.2*k[1][j] + (546.0/625.0)*k[2][j]
149 + (54.0/625.0)*k[3][j] - (378.0/625.0)*k[4][j];
150 this->sys->der_func(t,dq);
151 for (
int j = 0 ; j < this->sys->numVars(); j++) k[5][j] = step*dq[j];
154 for (
int j = 0; j < this->sys->numVars(); j++)
157 qq[j] += (1.0/24.0)*k[0][j] + (5.0/48.0)*k[3][j] +
158 (27.0/56.0)*k[4][j] + (125.0/336.0)*k[5][j];
161 fabs(k[0][j]/8.0+2.0*k[2][j]/3.0+k[3][j]/16.0-27.0*k[4][j]/56.0
162 -125.0*k[5][j]/336.0));
rk_45(ode_system< X > *sys, double err_tol, double h_max)
Definition: adevs_rk_45.h:70
Definition: adevs_rk_45.h:43
~rk_45()
Destructor.
Definition: adevs_rk_45.h:81
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
double integrate(double *q, double h_lim)
Definition: adevs_rk_45.h:97
Definition: adevs_hybrid.h:45
void advance(double *q, double h)
Definition: adevs_rk_45.h:90
Definition: adevs_hybrid.h:367