31 #ifndef _adevs_corrected_euler_h_
32 #define _adevs_corrected_euler_h_
34 #include "adevs_hybrid.h"
54 double integrate(
double* q,
double h_lim);
55 void advance(
double* q,
double h);
65 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 < 2; i++) k[i] =
new double[sys->
numVars()];
74 dq =
new double[sys->
numVars()];
75 qq =
new double[sys->
numVars()];
82 delete [] t;
delete [] qq;
delete [] dq;
83 for (
int i = 0; i < 2; i++)
delete [] k[i];
90 while ((dt = integrate(q,h)) < h) h -= dt;
97 double err = DBL_MAX, h = std::min(h_cur*1.1,std::min(h_max,h_lim));
100 for (
int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
104 if (err <= err_tol) {
105 if (h_lim >= h_cur) h_cur = h;
110 double h_guess = 0.8*err_tol*h/fabs(err);
111 if (h < h_guess) h *= 0.8;
116 for (
int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
120 template <
typename X>
125 this->sys->der_func(qq,dq);
126 for (j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
128 for (j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.5*k[0][j];
129 this->sys->der_func(t,dq);
130 for (j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
133 for (j = 0; j < this->sys->numVars(); j++) {
135 err = std::max(err,fabs(k[0][j]-k[1][j]));
Definition: adevs_corrected_euler.h:43
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
~corrected_euler()
Destructor.
Definition: adevs_corrected_euler.h:80
double integrate(double *q, double h_lim)
Definition: adevs_corrected_euler.h:94
corrected_euler(ode_system< X > *sys, double err_tol, double h_max)
Definition: adevs_corrected_euler.h:69
Definition: adevs_hybrid.h:45
void advance(double *q, double h)
Definition: adevs_corrected_euler.h:87
Definition: adevs_hybrid.h:367