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,
101 h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
104 for (
int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
108 if (err <= err_tol) {
109 if (h_cur <= h_lim) h_cur = h;
114 double h_guess = 0.8*pow(err_tol*pow(h,4.0)/fabs(err),0.25);
115 if (h < h_guess) h *= 0.8;
120 for (
int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
124 template <
typename X>
128 this->sys->der_func(qq,dq);
129 for (
int j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
131 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.5*k[0][j];
132 this->sys->der_func(t,dq);
133 for (
int j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
135 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.25*(k[0][j]+k[1][j]);
136 this->sys->der_func(t,dq);
137 for (
int j = 0; j < this->sys->numVars(); j++) k[2][j] = step*dq[j];
139 for (
int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] - k[1][j] + 2.0*k[2][j];
140 this->sys->der_func(t,dq);
141 for (
int j = 0; j < this->sys->numVars(); j++) k[3][j] = step*dq[j];
143 for (
int j = 0; j < this->sys->numVars(); j++)
144 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];
145 this->sys->der_func(t,dq);
146 for (
int j = 0; j < this->sys->numVars(); j++) k[4][j] = step*dq[j];
148 for (
int j = 0; j < this->sys->numVars(); j++)
149 t[j] = qq[j] + (28.0/625.0)*k[0][j] - 0.2*k[1][j] + (546.0/625.0)*k[2][j]
150 + (54.0/625.0)*k[3][j] - (378.0/625.0)*k[4][j];
151 this->sys->der_func(t,dq);
152 for (
int j = 0 ; j < this->sys->numVars(); j++) k[5][j] = step*dq[j];
155 for (
int j = 0; j < this->sys->numVars(); j++)
158 qq[j] += (1.0/24.0)*k[0][j] + (5.0/48.0)*k[3][j] +
159 (27.0/56.0)*k[4][j] + (125.0/336.0)*k[5][j];
162 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
163 -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
Definition: adevs_fmi.h:57
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:377