ode_rk8pd_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* ode-initval/rk8pd.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 #ifndef O2SCL_GSL_RK8PD_H
43 #define O2SCL_GSL_RK8PD_H
44 
45 /** \file ode_rk8pd_gsl.h
46  \brief File defining \ref o2scl::ode_rk8pd_gsl
47 */
48 
49 #include <o2scl/ode_step.h>
50 #include <o2scl/err_hnd.h>
51 
52 #ifndef DOXYGEN_NO_O2NS
53 namespace o2scl {
54 #endif
55 
56  /** \brief Embedded Runge-Kutta Prince-Dormand ODE stepper (GSL)
57 
58  Based on \ref Prince81 .
59 
60  There is an example for the usage of this class in
61  <tt>examples/ex_ode.cpp</tt> documented in the \ref ex_ode_sect
62  section.
63  */
64  template<class vec_y_t=boost::numeric::ublas::vector<double>,
65  class vec_dydx_t=vec_y_t, class vec_yerr_t=vec_y_t,
66  class func_t=ode_funct11>
67  class ode_rk8pd_gsl : public ode_step<vec_y_t,
68  vec_dydx_t,vec_yerr_t,func_t> {
69 
70  protected:
71 
72 #ifdef O2SCL_NEVER_DEFINED
73  }{
74 #endif
75 
76  /// \name Storage for the intermediate steps
77  //@{
78  vec_dydx_t k2, k3, k4, k5, k6, k7;
79  vec_dydx_t k8, k9, k10, k11, k12, k13;
80  vec_y_t ytmp;
81  //@}
82 
83  /// Size of allocated vectors
84  size_t ndim;
85 
86  /** \name Storage for the coefficients
87  */
88  //@{
89  double Abar[13], A[12], ah[10], b21, b3[2], b4[3], b5[4], b6[5];
90  double b7[6], b8[7], b9[8], b10[9], b11[10], b12[11], b13[12];
91  //@}
92 
93  public:
94 
95  ode_rk8pd_gsl() {
96  this->order=8;
97 
98  Abar[0]=14005451.0/335480064.0;
99  Abar[1]=0.0;
100  Abar[2]=0.0;
101  Abar[3]=0.0;
102  Abar[4]=0.0;
103  Abar[5]=-59238493.0/1068277825.0;
104  Abar[6]=181606767.0/758867731.0;
105  Abar[7]=561292985.0/797845732.0;
106  Abar[8]=-1041891430.0/1371343529.0;
107  Abar[9]=760417239.0/1151165299.0;
108  Abar[10]=118820643.0/751138087.0;
109  Abar[11]=-528747749.0/2220607170.0;
110  Abar[12]=1.0/4.0;
111 
112  A[0]=13451932.0/455176623.0;
113  A[1]=0.0;
114  A[2]=0.0;
115  A[3]=0.0;
116  A[4]=0.0;
117  A[5]=-808719846.0/976000145.0;
118  A[6]=1757004468.0/5645159321.0;
119  A[7]=656045339.0/265891186.0;
120  A[8]=-3867574721.0/1518517206.0;
121  A[9]=465885868.0/322736535.0;
122  A[10]=53011238.0/667516719.0;
123  A[11]=2.0/45.0;
124 
125  ah[0]=1.0/18.0;
126  ah[1]=1.0/12.0;
127  ah[2]=1.0/8.0;
128  ah[3]=5.0/16.0;
129  ah[4]=3.0/8.0;
130  ah[5]=59.0/400.0;
131  ah[6]=93.0/200.0;
132  ah[7]=5490023248.0/9719169821.0;
133  ah[8]=13.0/20.0;
134  ah[9]=1201146811.0/1299019798.0;
135 
136  b21=1.0/18.0;
137 
138  b3[0]=1.0/48.0;
139  b3[1]=1.0/16.0;
140 
141  b4[0]=1.0/32.0;
142  b4[1]=0.0;
143  b4[2]=3.0/32.0;
144 
145  b5[0]=5.0/16.0;
146  b5[1]=0.0;
147  b5[2]=-75.0/64.0;
148  b5[3]=75.0/64.0;
149 
150  b6[0]=3.0/80.0;
151  b6[1]=0.0;
152  b6[2]=0.0;
153  b6[3]=3.0/16.0;
154  b6[4]=3.0/20.0;
155 
156  b7[0]=29443841.0/614563906.0;
157  b7[1]=0.0;
158  b7[2]=0.0;
159  b7[3]=77736538.0/692538347.0;
160  b7[4]=-28693883.0/1125000000.0;
161  b7[5]=23124283.0/1800000000.0;
162 
163  b8[0]=16016141.0/946692911.0;
164  b8[1]=0.0;
165  b8[2]=0.0;
166  b8[3]=61564180.0/158732637.0;
167  b8[4]=22789713.0/633445777.0;
168  b8[5]=545815736.0/2771057229.0;
169  b8[6]=-180193667.0/1043307555.0;
170 
171  b9[0]=39632708.0/573591083.0;
172  b9[1]=0.0;
173  b9[2]=0.0;
174  b9[3]=-433636366.0/683701615.0;
175  b9[4]=-421739975.0/2616292301.0;
176  b9[5]=100302831.0/723423059.0;
177  b9[6]=790204164.0/839813087.0;
178  b9[7]=800635310.0/3783071287.0;
179 
180  b10[0]=246121993.0/1340847787.0;
181  b10[1]=0.0;
182  b10[2]=0.0;
183  b10[3]=-37695042795.0/15268766246.0;
184  b10[4]=-309121744.0/1061227803.0;
185  b10[5]=-12992083.0/490766935.0;
186  b10[6]=6005943493.0/2108947869.0;
187  b10[7]=393006217.0/1396673457.0;
188  b10[8]=123872331.0/1001029789.0;
189 
190  b11[0]=-1028468189.0/846180014.0;
191  b11[1]=0.0;
192  b11[2]=0.0;
193  b11[3]=8478235783.0/508512852.0;
194  b11[4]=1311729495.0/1432422823.0;
195  b11[5]=-10304129995.0/1701304382.0;
196  b11[6]=-48777925059.0/3047939560.0;
197  b11[7]=15336726248.0/1032824649.0;
198  b11[8]=-45442868181.0/3398467696.0;
199  b11[9]=3065993473.0/597172653.0;
200 
201  b12[0]=185892177.0/718116043.0;
202  b12[1]=0.0;
203  b12[2]=0.0;
204  b12[3]=-3185094517.0/667107341.0;
205  b12[4]=-477755414.0/1098053517.0;
206  b12[5]=-703635378.0/230739211.0;
207  b12[6]=5731566787.0/1027545527.0;
208  b12[7]=5232866602.0/850066563.0;
209  b12[8]=-4093664535.0/808688257.0;
210  b12[9]=3962137247.0/1805957418.0;
211  b12[10]=65686358.0/487910083.0;
212 
213  b13[0]=403863854.0/491063109.0;
214  b13[1]=0.0;
215  b13[2]=0.0;
216  b13[3]=-5068492393.0/434740067.0;
217  b13[4]=-411421997.0/543043805.0;
218  b13[5]=652783627.0/914296604.0;
219  b13[6]=11173962825.0/925320556.0;
220  b13[7]=-13158990841.0/6184727034.0;
221  b13[8]=3936647629.0/1978049680.0;
222  b13[9]=-160528059.0/685178525.0;
223  b13[10]=248638103.0/1413531060.0;
224  b13[11]=0.0;
225 
226  ndim=0;
227  }
228 
229  virtual ~ode_rk8pd_gsl() {
230  }
231 
232  /** \brief Perform an integration step
233 
234  Given initial value of the n-dimensional function in \c y and
235  the derivative in \c dydx (which must be computed beforehand)
236  at the point \c x, take a step of size \c h giving the result
237  in \c yout, the uncertainty in \c yerr, and the new derivative
238  in \c dydx_out using function \c derivs to calculate
239  derivatives. The parameters \c yout and \c y and the
240  parameters \c dydx_out and \c dydx may refer to the same
241  object.
242 
243  If \c derivs always returns zero, then this function will
244  also return zero. If not, <tt>step()</tt> will return the first
245  non-zero value which was obtained in a call to \c derivs .
246  The error handler is never called.
247  */
248  virtual int step(double x, double h, size_t n, vec_y_t &y,
249  vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr,
250  vec_dydx_t &dydx_out, func_t &derivs) {
251 
252  int ret=0;
253  size_t i;
254 
255  if (ndim!=n) {
256  k2.resize(n);
257  k3.resize(n);
258  k4.resize(n);
259  k5.resize(n);
260  k6.resize(n);
261  k7.resize(n);
262  k8.resize(n);
263  k9.resize(n);
264  k10.resize(n);
265  k11.resize(n);
266  k12.resize(n);
267  k13.resize(n);
268  ytmp.resize(n);
269 
270  ndim=n;
271  }
272 
273  for (i=0;i<n;i++) {
274  ytmp[i]=y[i]+b21*h*dydx[i];
275  }
276 
277  error_update(ret,derivs(x+ah[0]*h,n,ytmp,k2));
278 
279  for (i=0;i<n;i++) {
280  ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
281  }
282 
283  error_update(ret,derivs(x+ah[1]*h,n,ytmp,k3));
284 
285  for (i=0;i<n;i++) {
286  ytmp[i]=y[i]+h*(b4[0]*dydx[i]+b4[2]*k3[i]);
287  }
288 
289  error_update(ret,derivs(x+ah[2]*h,n,ytmp,k4));
290 
291  for (i=0;i<n;i++) {
292  ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[2]*k3[i]+b5[3]*k4[i]);
293  }
294 
295  error_update(ret,derivs(x+ah[3]*h,n,ytmp,k5));
296 
297  for (i=0;i<n;i++) {
298  ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[3]*k4[i]+b6[4]*k5[i]);
299  }
300 
301  error_update(ret,derivs(x+ah[4]*h,n,ytmp,k6));
302 
303  for (i=0;i<n;i++) {
304  ytmp[i]=y[i]+h*(b7[0]*dydx[i]+b7[3]*k4[i]+b7[4]*k5[i]+b7[5]*k6[i]);
305  }
306 
307  error_update(ret,derivs(x+ah[5]*h,n,ytmp,k7));
308 
309  for (i=0;i<n;i++) {
310  ytmp[i]=y[i]+h*(b8[0]*dydx[i]+b8[3]*k4[i]+b8[4]*k5[i]+b8[5]*k6[i]+
311  b8[6]*k7[i]);
312  }
313 
314  error_update(ret,derivs(x+ah[6]*h,n,ytmp,k8));
315 
316  for (i=0;i<n;i++) {
317  ytmp[i]=y[i]+h*(b9[0]*dydx[i]+b9[3]*k4[i]+b9[4]*k5[i]+b9[5]*k6[i]+
318  b9[6]*k7[i]+b9[7]*k8[i]);
319  }
320 
321  error_update(ret,derivs(x+ah[7]*h,n,ytmp,k9));
322 
323  for (i=0;i<n;i++) {
324  ytmp[i]=y[i]+h*(b10[0]*dydx[i]+b10[3]*k4[i]+b10[4]*k5[i]+
325  b10[5]*k6[i]+b10[6]*k7[i]+b10[7]*k8[i]+
326  b10[8]*k9[i]);
327  }
328 
329  error_update(ret,derivs(x+ah[8]*h,n,ytmp,k10));
330 
331  for (i=0;i<n;i++) {
332  ytmp[i]=y[i]+h*(b11[0]*dydx[i]+b11[3]*k4[i]+b11[4]*k5[i]+
333  b11[5]*k6[i]+b11[6]*k7[i]+b11[7]*k8[i]+
334  b11[8]*k9[i]+b11[9]*k10[i]);
335  }
336 
337  error_update(ret,derivs(x+ah[9]*h,n,ytmp,k11));
338 
339  for (i=0;i<n;i++) {
340  ytmp[i]=y[i]+h*(b12[0]*dydx[i]+b12[3]*k4[i]+b12[4]*k5[i]+
341  b12[5]*k6[i]+b12[6]*k7[i]+b12[7]*k8[i]+
342  b12[8]*k9[i]+b12[9]*k10[i]+b12[10]*k11[i]);
343  }
344 
345  error_update(ret,derivs(x+h,n,ytmp,k12));
346 
347  for (i=0;i<n;i++) {
348  ytmp[i]=y[i]+h*(b13[0]*dydx[i]+b13[3]*k4[i]+b13[4]*k5[i]+
349  b13[5]*k6[i]+b13[6]*k7[i]+b13[7]*k8[i]+
350  b13[8]*k9[i]+b13[9]*k10[i]+b13[10]*k11[i]+
351  b13[11]*k12[i]);
352  }
353 
354  error_update(ret,derivs(x+h,n,ytmp,k13));
355 
356  // final sum
357 
358  for (i=0;i<n;i++) {
359  double ksum8=Abar[0]*dydx[i]+Abar[5]*k6[i]+Abar[6]*k7[i]+
360  Abar[7]*k8[i]+Abar[8]*k9[i]+Abar[9]*k10[i]+
361  Abar[10]*k11[i]+Abar[11]*k12[i]+Abar[12]*k13[i];
362 
363  yout[i]=y[i]+h*ksum8;
364  }
365 
366  // We put this before the last function evaluation, in contrast
367  // to the GSL version, so that the dydx[i] that appears in the
368  // for loop below isn't modified by the subsequent derivative
369  // evaluation using dydx_out. (The user could have given the
370  // same vector for both)
371  for (i=0;i<n;i++) {
372 
373  double ksum8=Abar[0]*dydx[i]+Abar[5]*k6[i]+Abar[6]*k7[i]+
374  Abar[7]*k8[i]+Abar[8]*k9[i]+Abar[9]*k10[i]+
375  Abar[10]*k11[i]+Abar[11]*k12[i]+Abar[12]*k13[i];
376  double ksum7=A[0]*dydx[i]+A[5]*k6[i]+A[6]*k7[i]+A[7]*k8[i]+
377  A[8]*k9[i]+A[9]*k10[i]+A[10]*k11[i]+A[11]*k12[i];
378 
379  yerr[i]=h*(ksum7-ksum8);
380  }
381 
382  error_update(ret,derivs(x+h,n,yout,dydx_out));
383 
384  return ret;
385  }
386 
387  };
388 
389 #ifndef DOXYGEN_NO_O2NS
390 }
391 #endif
392 
393 #endif
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
Embedded Runge-Kutta Prince-Dormand ODE stepper (GSL)
Definition: ode_rk8pd_gsl.h:67
void error_update(int &ret, int err)
Update an error value err with the value in ret.
Definition: err_hnd.h:319
int order
The order of the ODE stepper.
Definition: ode_step.h:49
virtual int step(double x, double h, size_t n, vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out, func_t &derivs)
Perform an integration step.
size_t ndim
Size of allocated vectors.
Definition: ode_rk8pd_gsl.h:84
ODE stepper base [abstract base].
Definition: ode_step.h:41

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).