31 #ifndef _adevs_corrected_euler_h_ 32 #define _adevs_corrected_euler_h_ 35 #include "adevs_hybrid.h" 55 double integrate(
double* q,
double h_lim);
56 void advance(
double* q,
double h);
66 double trial_step(
double h);
72 ode_solver<X>(sys),err_tol(err_tol),h_max(h_max),h_cur(h_max)
74 for (
int i = 0; i < 2; i++) k[i] =
new double[sys->
numVars()];
75 dq =
new double[sys->
numVars()];
76 qq =
new double[sys->
numVars()];
83 delete [] t;
delete [] qq;
delete [] dq;
84 for (
int i = 0; i < 2; i++)
delete [] k[i];
91 while ((dt =
integrate(q,h)) < h) h -= dt;
99 h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
102 for (
int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
106 if (err <= err_tol) {
107 if (h_lim >= h_cur) h_cur = h;
112 double h_guess = 0.8*err_tol*h/fabs(err);
113 if (h < h_guess) h *= 0.8;
118 for (
int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
122 template <
typename X>
127 this->sys->der_func(qq,dq);
128 for (j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
130 for (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 (j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
135 for (j = 0; j < this->sys->numVars(); j++) {
137 err = std::max<double>(err,fabs(k[0][j]-k[1][j]));
Definition: adevs_corrected_euler.h:44
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
Definition: adevs_fmi.h:57
~corrected_euler()
Destructor.
Definition: adevs_corrected_euler.h:81
double integrate(double *q, double h_lim)
Definition: adevs_corrected_euler.h:95
corrected_euler(ode_system< X > *sys, double err_tol, double h_max)
Definition: adevs_corrected_euler.h:70
Definition: adevs_hybrid.h:45
void advance(double *q, double h)
Definition: adevs_corrected_euler.h:88
Definition: adevs_hybrid.h:377