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