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

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