adevs
adevs_corrected_euler.h
1 
31 #ifndef _adevs_corrected_euler_h_
32 #define _adevs_corrected_euler_h_
33 #include <cmath>
34 #include <algorithm>
35 #include "adevs_hybrid.h"
36 
37 namespace adevs
38 {
39 
44 template <typename X> class corrected_euler:
45  public ode_solver<X>
46 {
47  public:
52  corrected_euler(ode_system<X>* sys, double err_tol, double h_max);
55  double integrate(double* q, double h_lim);
56  void advance(double* q, double h);
57  private:
58  double *dq, // derivative
59  *qq, // trial solution
60  *t, // temporary variable for computing k2
61  *k[2]; // k1 and k2
62  const double err_tol; // Error tolerance
63  const double h_max; // Maximum time step
64  double h_cur; // Previous time step that satisfied error constraint
65  // Compute a step of size h, put it in qq, and return the error
66  double trial_step(double h);
67 };
68 
69 template <typename X>
71  double h_max):
72  ode_solver<X>(sys),err_tol(err_tol),h_max(h_max),h_cur(h_max)
73 {
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()];
77  t = new double[sys->numVars()];
78 }
79 
80 template <typename X>
82 {
83  delete [] t; delete [] qq; delete [] dq;
84  for (int i = 0; i < 2; i++) delete [] k[i];
85 }
86 
87 template <typename X>
88 void corrected_euler<X>::advance(double* q, double h)
89 {
90  double dt;
91  while ((dt = integrate(q,h)) < h) h -= dt;
92 }
93 
94 template <typename X>
95 double corrected_euler<X>::integrate(double* q, double h_lim)
96 {
97  // Initial error estimate and step size
98  double err = DBL_MAX,
99  h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
100  for (;;) {
101  // Copy q to the trial vector
102  for (int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
103  // Make the trial step which will be stored in qq
104  err = trial_step(h);
105  // If the error is ok, then we have found the proper step size
106  if (err <= err_tol) { // Keep h if shrunk to control the error
107  if (h_lim >= h_cur) h_cur = h;
108  break;
109  }
110  // Otherwise shrink the step size and try again
111  else {
112  double h_guess = 0.8*err_tol*h/fabs(err);
113  if (h < h_guess) h *= 0.8;
114  else h = h_guess;
115  }
116  }
117  // Put the trial solution in q and return the selected step size
118  for (int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
119  return h;
120 }
121 
122 template <typename X>
123 double corrected_euler<X>::trial_step(double step)
124 {
125  int j;
126  // Compute k1
127  this->sys->der_func(qq,dq);
128  for (j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
129  // Compute k2
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];
133  // Compute next state and approximate error
134  double err = 0.0;
135  for (j = 0; j < this->sys->numVars(); j++) {
136  qq[j] += k[1][j]; // Next state
137  err = std::max<double>(err,fabs(k[0][j]-k[1][j])); // Maximum error
138  }
139  return err; // Return the error
140 }
141 
142 } // end of namespace
143 #endif
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