adevs
adevs_rk_45.h
1 
31 #ifndef _adevs_rk_45_h_
32 #define _adevs_rk_45_h_
33 #include "adevs_hybrid.h"
34 #include <cmath>
35 
36 namespace adevs
37 {
38 
43 template <typename X> class rk_45:
44  public ode_solver<X>
45 {
46  public:
52  rk_45(ode_system<X>* sys, double err_tol, double h_max);
54  ~rk_45();
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 variables for computing stages
61  *k[6]; // the six RK stages
62  const double err_tol; // Error tolerance
63  const double h_max; // Maximum time step
64  double h_cur; // Previous successful step size
65  // Compute a trial step of size h, store the result in qq, and return the error
66  double trial_step(double h);
67 };
68 
69 template <typename X>
70 rk_45<X>::rk_45(ode_system<X>* sys, double err_tol, 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 < 6; i++)
74  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 [] dq;
84  delete [] t;
85  for (int i = 0; i < 6; i++)
86  delete [] k[i];
87 }
88 
89 template <typename X>
90 void rk_45<X>::advance(double* q, double h)
91 {
92  double dt;
93  while ((dt = integrate(q,h)) < h) h -= dt;
94 }
95 
96 template <typename X>
97 double rk_45<X>::integrate(double* q, double h_lim)
98 {
99  // Initial error estimate and step size
100  double err = DBL_MAX,
101  h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
102  for (;;) {
103  // Copy q to the trial vector
104  for (int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
105  // Make the trial step which will be stored in qq
106  err = trial_step(h);
107  // If the error is ok, then we have found the proper step size
108  if (err <= err_tol) {
109  if (h_cur <= h_lim) h_cur = h;
110  break;
111  }
112  // Otherwise shrink the step size and try again
113  else {
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;
116  else h = h_guess;
117  }
118  }
119  // Copy the trial solution to q and return the step size that was selected
120  for (int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
121  return h;
122 }
123 
124 template <typename X>
125 double rk_45<X>::trial_step(double step)
126 {
127  // Compute k1
128  this->sys->der_func(qq,dq);
129  for (int j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
130  // Compute k2
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];
134  // Compute k3
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];
138  // Compute k4
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];
142  // Compute k5
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];
147  // Compute k6
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];
153  // Compute next state and the approximate error
154  double err = 0.0;
155  for (int j = 0; j < this->sys->numVars(); j++)
156  {
157  // Next state
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];
160  // Componennt wise maximum of the approximate error
161  err = std::max(err,
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));
164  }
165  // Return the error
166  return err;
167 }
168 
169 } // end of namespace
170 #endif
171 
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