ode_rkck_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/rkck.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_RKCK_H
43 #define O2SCL_GSL_RKCK_H
44 
45 /** \file ode_rkck_gsl.h
46  \brief File defining \ref o2scl::ode_rkck_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 Cash-Karp embedded Runge-Kutta ODE stepper (GSL)
58 
59  Based on \ref Cash90 .
60 
61  There is an example for the usage of this class in
62  <tt>examples/ex_ode.cpp</tt> documented in the \ref ex_ode_sect
63  section.
64  */
65  template<class vec_y_t=boost::numeric::ublas::vector<double>,
66  class vec_dydx_t=vec_y_t, class vec_yerr_t=vec_y_t,
67  class func_t=ode_funct11>
68  class ode_rkck_gsl : public ode_step<vec_y_t,
69  vec_dydx_t,vec_yerr_t,func_t> {
70 
71  protected:
72 
73  /// \name Storage for the intermediate steps
74  //@{
75  vec_y_t ytmp;
76  vec_dydx_t k2, k3, k4, k5, k6;
77  //@}
78 
79  /// Size of allocated vectors
80  size_t ndim;
81 
82  /** \name Storage for the coefficients
83  */
84  //@{
85  double ah[5], b3[2], b4[3], b5[4], b6[5], ec[7];
86  double b21, c1, c3, c4, c6;
87  //@}
88 
89  public:
90 
91  ode_rkck_gsl() {
92  this->order=5;
93 
94  ah[0]=1.0/5.0;
95  ah[1]=3.0/10.0;
96  ah[2]=3.0/5.0;
97  ah[3]=1.0;
98  ah[4]=7.0/8.0;
99 
100  b3[0]=3.0/40.0;
101  b3[1]=9.0/40.0;
102 
103  b4[0]=3.0/10.0;
104  b4[1]=-9.0/10.0;
105  b4[2]=12.0/10.0;
106 
107  b5[0]=-11.0/54.0;
108  b5[1]=5.0/2.0;
109  b5[2]=-70.0/27.0;
110  b5[3]=35.0/27.0;
111 
112  b6[0]=1631.0/55296.0;
113  b6[1]=175.0/512.0;
114  b6[2]=575.0/13824.0;
115  b6[3]=44275.0/110592.0;
116  b6[4]=253.0/4096.0;
117 
118  ec[0]=0.0;
119  ec[1]=37.0/378.0-2825.0/27648.0;
120  ec[2]=0.0;
121  ec[3]=250.0/621.0-18575.0/48384.0;
122  ec[4]=125.0/594.0-13525.0/55296.0;
123  ec[5]=-277.0/14336.0;
124  ec[6]=512.0/1771.0-1.0/4.0;
125 
126  b21=1.0/5.0;
127 
128  c1=37.0/378.0;
129  c3=250.0/621.0;
130  c4=125.0/594.0;
131  c6=512.0/1771.0;
132 
133  ndim=0;
134  }
135 
136  virtual ~ode_rkck_gsl() {
137  }
138 
139  /** \brief Perform an integration step
140 
141  Given initial value of the n-dimensional function in \c y and
142  the derivative in \c dydx (which must be computed beforehand)
143  at the point \c x, take a step of size \c h giving the result
144  in \c yout, the uncertainty in \c yerr, and the new derivative
145  in \c dydx_out using function \c derivs to calculate
146  derivatives. The parameters \c yout and \c y and the
147  parameters \c dydx_out and \c dydx may refer to the same
148  object.
149 
150  If \c derivs always returns zero, then this function will
151  also return zero. If not, <tt>step()</tt> will return the first
152  non-zero value which was obtained in a call to \c derivs .
153  The error handler is never called.
154  */
155  virtual int step(double x, double h, size_t n, vec_y_t &y,
156  vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr,
157  vec_dydx_t &dydx_out, func_t &derivs) {
158 
159  int ret=0;
160  size_t i;
161 
162  if (ndim!=n) {
163  k2.resize(n);
164  k3.resize(n);
165  k4.resize(n);
166  k5.resize(n);
167  k6.resize(n);
168  ytmp.resize(n);
169 
170  ndim=n;
171  }
172 
173  for (i=0;i<n;i++) {
174  ytmp[i]=y[i]+b21*h*dydx[i];
175  }
176 
177  o2scl::error_update(ret,derivs(x+ah[0]*h,n,ytmp,k2));
178 
179  for (i=0;i<n;i++) {
180  ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
181  }
182 
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  o2scl::error_update(ret,derivs(x+ah[2]*h,n,ytmp,k4));
190 
191  for (i=0;i<n;i++) {
192  ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[1]*k2[i]+b5[2]*k3[i]+
193  b5[3]*k4[i]);
194  }
195 
196  o2scl::error_update(ret,derivs(x+ah[3]*h,n,ytmp,k5));
197 
198  for (i=0;i<n;i++) {
199  ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[1]*k2[i]+b6[2]*k3[i]+
200  b6[3]*k4[i]+b6[4]*k5[i]);
201  }
202 
203  o2scl::error_update(ret,derivs(x+ah[4]*h,n,ytmp,k6));
204 
205  for (i=0;i<n;i++) {
206  yout[i]=y[i]+h*(c1*dydx[i]+c3*k3[i]+c4*k4[i]+c6*k6[i]);
207  }
208 
209  // We put this before the last function evaluation, in contrast
210  // to the GSL version, so that the dydx[i] that appears in the
211  // for loop below isn't modified by the subsequent derivative
212  // evaluation using dydx_out. (The user could have given the
213  // same vector for both)
214  for (i=0;i<n;i++) {
215  yerr[i]=h*(ec[1]*dydx[i]+ec[3]*k3[i]+ec[4]*k4[i]+ec[5]*k5[i]+
216  ec[6]*k6[i]);
217  }
218 
219  o2scl::error_update(ret,derivs(x+h,n,yout,dydx_out));
220 
221  return ret;
222  }
223 
224  };
225 
226 #ifndef DOXYGEN_NO_O2NS
227 }
228 #endif
229 
230 #endif
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.
Definition: ode_rkck_gsl.h:155
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
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
Cash-Karp embedded Runge-Kutta ODE stepper (GSL)
Definition: ode_rkck_gsl.h:68
size_t ndim
Size of allocated vectors.
Definition: ode_rkck_gsl.h:80
ODE stepper base [abstract base].
Definition: ode_step.h:41

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