ode_rkf45_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/rkf45.c
24  *
25  * Copyright (C) 2001, 2004, 2007 Brian Gough
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_RKF45_H
43 #define O2SCL_GSL_RKF45_H
44 
45 /** \file ode_rkf45_gsl.h
46  \brief File defining \ref o2scl::ode_rkf45_gsl
47 */
48 
49 #include <o2scl/err_hnd.h>
50 #include <o2scl/ode_funct.h>
51 #include <o2scl/ode_step.h>
52 
53 #ifndef DOXYGEN_NO_O2NS
54 namespace o2scl {
55 #endif
56 
57  /** \brief Runge-Kutta-Fehlberg embedded Runge-Kutta ODE stepper (GSL)
58 
59  Based on \ref Hairer09 .
60 
61  \todo Check this because it may not give exact dydt_out.
62  */
63  template<class vec_y_t=boost::numeric::ublas::vector<double>,
64  class vec_dydx_t=vec_y_t, class vec_yerr_t=vec_y_t,
65  class func_t=ode_funct11>
66  class ode_rkf45_gsl : public ode_step<vec_y_t,
67  vec_dydx_t,vec_yerr_t,func_t> {
68 
69  protected:
70 
71  /// \name Storage for the intermediate steps
72  //@{
73  vec_dydx_t k2, k3, k4, k5, k6;
74  vec_y_t ytmp;
75  //@}
76 
77  /// Size of allocated vectors
78  size_t ndim;
79 
80  /** \name Storage for the coefficients
81  */
82  //@{
83  double ah[5], b3[2], b4[3], b5[4], b6[5];
84  double c1, c3, c4, c5, c6;
85  double ec[7];
86  //@}
87 
88  public:
89 
90  ode_rkf45_gsl() {
91  this->order=5;
92 
93  ah[0]=1.0/4.0;
94  ah[1]=3.0/8.0;
95  ah[2]=12.0/13.0;
96  ah[3]=1.0;
97  ah[4]=1.0/2.0;
98 
99  b3[0]=3.0/32.0;
100  b3[1]=9.0/32.0;
101 
102  b4[0]=1932.0/2197.0;
103  b4[1]=-7200.0/2197.0;
104  b4[2]=7296.0/2197.0;
105 
106  b5[0]=8341.0/4104.0;
107  b5[1]=-32832.0/4104.0;
108  b5[2]=29440.0/4104.0;
109  b5[3]=-845.0/4104.0;
110 
111  b6[0]=-6080.0/20520.0;
112  b6[1]=41040.0/20520.0;
113  b6[2]=-28352.0/20520.0;
114  b6[3]=9295.0/20520.0;
115  b6[4]=-5643.0/20520.0;
116 
117  c1=902880.0/7618050.0;
118  c3=3953664.0/7618050.0;
119  c4=3855735.0/7618050.0;
120  c5=-1371249.0/7618050.0;
121  c6=277020.0/7618050.0;
122 
123  ec[0]=0.0;
124  ec[1]=1.0/360.0;
125  ec[2]=0.0;
126  ec[3]=-128.0/4275.0;
127  ec[4]=-2197.0/75240.0;
128  ec[5]=1.0/50.0;
129  ec[6]=2.0/55.0;
130 
131  ndim=0;
132  }
133 
134  virtual ~ode_rkf45_gsl() {
135  }
136 
137  /** \brief Perform an integration step
138 
139  Given initial value of the n-dimensional function in \c y and
140  the derivative in \c dydx (which must be computed beforehand) at
141  the point \c x, take a step of size \c h giving the result in \c
142  yout, the uncertainty in \c yerr, and the new derivative in \c
143  dydx_out using function \c derivs to calculate derivatives. The
144  parameters \c yout and \c y and the parameters \c dydx_out and
145  \c dydx may refer to the same object.
146 
147  If \c derivs always returns zero, then this function will
148  also return zero. If not, <tt>step()</tt> will return the first
149  non-zero value which was obtained in a call to \c derivs .
150  The error handler is never called.
151  */
152  virtual int step(double x, double h, size_t n, vec_y_t &y, vec_dydx_t &dydx,
153  vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out,
154  func_t &derivs) {
155 
156  int ret=0;
157  size_t i;
158 
159  if (ndim!=n) {
160  k2.resize(n);
161  k3.resize(n);
162  k4.resize(n);
163  k5.resize(n);
164  k6.resize(n);
165  ytmp.resize(n);
166 
167  ndim=n;
168  }
169 
170  // k1 step
171  for (i=0;i<n;i++) {
172  ytmp[i]=y[i]+ah[0]*h*dydx[i];
173  }
174 
175  // k2 step
176  o2scl::error_update(ret,derivs(x+ah[0]*h,n,ytmp,k2));
177 
178  for (i=0;i<n;i++) {
179  ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
180  }
181 
182  // k3 step
183  o2scl::error_update(ret,derivs(x+ah[1]*h,n,ytmp,k3));
184 
185  for (i=0;i<n;i++) {
186  ytmp[i]=y[i]+h*(b4[0]*dydx[i]+b4[1]*k2[i]+b4[2]*k3[i]);
187  }
188 
189  // k4 step
190  o2scl::error_update(ret,derivs(x+ah[2]*h,n,ytmp,k4));
191 
192  for (i=0;i<n;i++) {
193  ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[1]*k2[i]+b5[2]*k3[i]+
194  b5[3]*k4[i]);
195  }
196 
197  // k5 step
198  o2scl::error_update(ret,derivs(x+ah[3]*h,n,ytmp,k5));
199 
200  for (i=0;i<n;i++) {
201  ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[1]*k2[i]+b6[2]*k3[i]+
202  b6[3]*k4[i]+b6[4]*k5[i]);
203  }
204 
205  // k6 step and final sum
206  o2scl::error_update(ret,derivs(x+ah[4]*h,n,ytmp,k6));
207 
208  for (i=0;i<n;i++) {
209  yout[i]=y[i]+h*(c1*dydx[i]+c3*k3[i]+c4*k4[i]+c5*k5[i]+c6*k6[i]);
210  }
211 
212  // We put this before the last function evaluation, in contrast
213  // to the GSL version, so that the dydx[i] that appears in the
214  // for loop below isn't modified by the subsequent derivative
215  // evaluation using dydx_out. (The user could have given the
216  // same vector for both)
217  for (i=0;i<n;i++) {
218  yerr[i]=h*(ec[1]*dydx[i]+ec[3]*k3[i]+ec[4]*k4[i]+ec[5]*k5[i]+
219  ec[6]*k6[i]);
220  }
221 
222  o2scl::error_update(ret,derivs(x+h,n,yout,dydx_out));
223 
224  return ret;
225  }
226 
227  };
228 
229 #ifndef DOXYGEN_NO_O2NS
230 }
231 #endif
232 
233 #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
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.
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
size_t ndim
Size of allocated vectors.
Definition: ode_rkf45_gsl.h:78
Runge-Kutta-Fehlberg embedded Runge-Kutta ODE stepper (GSL)
Definition: ode_rkf45_gsl.h:66
ODE stepper base [abstract base].
Definition: ode_step.h:41

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