astep_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-initval2/evolve.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_ASTEP_GSL_H
43 #define O2SCL_ASTEP_GSL_H
44 
45 /** \file astep_gsl.h
46  \brief File defining \ref o2scl::astep_gsl
47 */
48 
49 #include <string>
50 
51 #include <gsl/gsl_math.h>
52 #include <gsl/gsl_odeiv.h>
53 
54 #include <o2scl/astep.h>
55 #include <o2scl/vector.h>
56 #include <o2scl/string_conv.h>
57 
58 #ifndef DOXYGEN_NO_O2NS
59 namespace o2scl {
60 #endif
61 
62  /** \brief Control structure for astep_gsl
63 
64  This class implements both the "standard" and "scaled" step
65  control methods from GSL. The standard control method is the
66  default. To use the scaled control, set \ref standard to
67  <tt>false</tt> and set the scale for each component using
68  set_scale().
69 
70  The control object is a four parameter heuristic based on
71  absolute and relative errors \ref eps_abs and \ref eps_rel, and
72  scaling factors \ref a_y and \ref a_dydt for the system state
73  \f$ y(t) \f$ and derivatives \f$ y^{\prime}(t) \f$ respectively.
74 
75  The step-size adjustment procedure for this method begins by
76  computing the desired error level \f$ D_i \f$ for each
77  component. In the unscaled version,
78  \f[
79  D_i = \mathrm{eps\_abs}+\mathrm{eps\_rel} \times
80  \left( \mathrm{a\_y} | y_i| + \mathrm{a\_dydt}~h
81  | y_i^{\prime}| \right)
82  \f]
83  while in the scaled version the user specifies the scale for
84  each component, \f$ s_i \f$,
85  \f[
86  D_i = \mathrm{eps\_abs}~s_i+\mathrm{eps\_rel} \times
87  \left( \mathrm{a\_y} | y_i| + \mathrm{a\_dydt}~h
88  | y_i^{\prime}| \right)
89  \f]
90 
91  The desired error level \f$ D_i \f$ is compared to then observed
92  error \f$ E_i = |\mathrm{yerr}_i| \f$. If the observed error \f$
93  E \f$ exceeds the desired error level \f$ D \f$ by more than 10
94  percent for any component then the method reduces the step-size
95  by an appropriate factor,
96  \f[
97  h_{\mathrm{new}} = S~h_{\mathrm{old}}
98  \left(\frac{E}{D}\right)^{-1/q}
99  \f]
100  where \f$ q \f$ is the consistency order of the method (e.g. \f$
101  q=4 \f$ for 4(5) embedded RK), and \f$ S \f$ is a safety factor
102  of 0.9. The ratio \f$ E/D \f$ is taken to be the maximum of the
103  ratios \f$ E_i/D_i \f$.
104 
105  If the observed error E is less than 50 percent of the desired
106  error level \f$ D \f$ for the maximum ratio \f$ E_i/D_i \f$ then
107  the algorithm takes the opportunity to increase the step-size to
108  bring the error in line with the desired level,
109  \f[
110  h_{\mathrm{new}} = S~h_{\mathrm{old}}
111  \left(\frac{E}{D}\right)^{-1/(q+1)}
112  \f]
113  This encompasses all the standard error scaling methods. To
114  avoid uncontrolled changes in the stepsize, the overall scaling
115  factor is limited to the range 1/5 to 5.
116 
117  If the user specified fewer scaling parameters than the number
118  of ODEs, then the scaling parameters are reused as follows. If
119  there are \f$ N \f$ ODEs and \f$ M \f$ scaling parameters, then
120  for \f$ i>M \f$, the ith scaling parameter \f$ s_i \f$ is set to
121  be \f$ s_{i\%M} \f$ . If the user selects the scaled control by
122  setting \ref standard to <tt>false</tt> and no scale parameters
123  are specified, this class reverts to the standard control.
124 
125  \todo Double check that the improvements in the ode-initval2
126  routines are available here
127  */
128  template<class vec_y_t=boost::numeric::ublas::vector<double>,
129  class vec_dydx_t=vec_y_t, class vec_yerr_t=vec_y_t>
131 
132  public:
133 
135 
136  protected:
137 
138  /// Scalings
139  ubvector scale_abs;
140 
141  public:
142 
143  /// \name Adjustment specification
144  //@{
145  /// No adjustment required
146  static const size_t hadj_nil=0;
147  /// Recommend step decrease
148  static const size_t hadj_dec=1;
149  /// Recommend step increase
150  static const size_t hadj_inc=2;
151  //@}
152 
153  /// Absolute precision (default \f$ 10^{-6} \f$)
154  double eps_abs;
155  /// Relative precision (default 0)
156  double eps_rel;
157  /// Function scaling factor (default 1)
158  double a_y;
159  /// Derivative scaling factor (default 0)
160  double a_dydt;
161  /// Use standard or scaled algorithm (default true)
162  bool standard;
163 
164  ode_control_gsl() {
165  eps_abs=1.0e-6;
166  eps_rel=0.0;
167  a_y=1.0;
168  a_dydt=0.0;
169  standard=true;
170  }
171 
172  virtual ~ode_control_gsl() {
173  }
174 
175  /** \brief Set the scaling for each differential equation
176  */
177  template<class svec_t> int set_scale(size_t nscal, const svec_t &scale) {
178  scale_abs.resize(nscal);
179  for(size_t i=0;i<nscal;i++) {
180  scale_abs[i]=scale[i];
181  }
182  return 0;
183  }
184 
185  virtual int hadjust(size_t dim, unsigned int ord, const vec_y_t &y,
186  vec_yerr_t &yerr, vec_dydx_t &yp, double &h) {
187 
188  const double S=0.9;
189  const double h_old=h;
190 
191  double rmax=DBL_MIN;
192  size_t i;
193 
194  for(i=0;i<dim;i++) {
195  double D0;
196  if (standard || scale_abs.size()==0) {
197  D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+eps_abs;
198  } else {
199  D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+
200  eps_abs*scale_abs[i%scale_abs.size()];
201  }
202  const double r=fabs(yerr[i]) / fabs(D0);
203  rmax=GSL_MAX_DBL(r, rmax);
204  }
205 
206  if (rmax > 1.1) {
207 
208  /* decrease step, no more than factor of 5, but a fraction S more
209  than scaling suggests (for better accuracy) */
210  double r= S / pow(rmax, 1.0/ord);
211 
212  if (r < 0.2) {
213  r=0.2;
214  }
215  h=r*h_old;
216  return hadj_dec;
217 
218  } else if (rmax < 0.5) {
219 
220  /* increase step, no more than factor of 5 */
221  double r=S / pow(rmax, 1.0/(ord+1.0));
222  if (r > 5.0) {
223  r=5.0;
224  }
225  /* don't allow any decrease caused by S<1 */
226  if (r < 1.0) {
227  r=1.0;
228  }
229  h=r*h_old;
230  return hadj_inc;
231 
232  } else {
233 
234  /* no change */
235  return hadj_nil;
236 
237  }
238 
239  }
240 
241  };
242 
243  /** \brief Adaptive ODE stepper (GSL)
244 
245  This class performs an adaptive step of a system of ODEs. To
246  modify the ODE stepper which is used, use the function \ref
247  astep_base::set_step().
248 
249  Note, this has been updated to correspond to the
250  <tt>ode-initval2</tt> functions in GSL.
251 
252  There is an example for the usage of this class in
253  <tt>examples/ex_ode.cpp</tt> documented in the \ref ex_ode_sect
254  section.
255 
256  \todo Document what happens when the stepper function returns
257  a non-zero value, as it's different now with the ode-initval2
258  function.
259  \todo Document count, failed_steps, etc.
260  \future Compare more directly to GSL
261 
262  Default template arguments
263  - \c func_t - \ref ode_funct11
264  - \c vec_t - \ref boost::numeric::ublas::vector < double >
265  */
266  template<class vec_y_t=boost::numeric::ublas::vector<double>,
267  class vec_dydx_t=vec_y_t, class vec_yerr_t=vec_y_t,
268  class func_t=ode_funct11 > class astep_gsl :
269  public astep_base<vec_y_t,vec_dydx_t,vec_yerr_t,func_t> {
270 
271  public:
272 
274 
275 #ifndef DOXYGEN_INTERNAL
276 
277  protected:
278 
279  /// Temporary storage for yout
280  vec_y_t yout_int;
281 
282  /// Internal storage for dydx
283  vec_dydx_t dydx_int;
284 
285  /// The size of the last step
286  double last_step;
287 
288  /// The number of steps
289  unsigned long int count;
290 
291  /// The number of failed steps
292  unsigned long int failed_steps;
293 
294  /// The size of the allocated vectors
295  size_t msize;
296 
297  /** \brief Apply the evolution for the next adaptive step
298 
299  This function is based on <tt>gsl_odeiv2_evolve_apply</tt>.
300 
301  \note This function requres that \c y, \c yout, \c dydx and \c
302  dydx_out are all distinct vectors.
303  */
304  int evolve_apply(double t0, double t1, double &t, double &h, size_t nvar,
305  vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout,
306  vec_yerr_t &yerr, vec_dydx_t &dydx_out,
307  func_t &derivs) {
308 
309  double h0=h;
310  int step_status;
311  bool final_step=false;
312  // Remaining step size, possibly less than h
313  double dt=t1-t0;
314 
315  if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0)) {
316  std::string str="Interval direction (t1-t0="+o2scl::dtos(dt)+
317  ") does not match step direction (h="+o2scl::dtos(h)+
318  " in astep_gsl::evolve_apply().";
319  O2SCL_ERR(str.c_str(),exc_einval);
320  }
321 
322  // Copy the initial point from y to yout. The stepper won't modify
323  // y, so we can undo the step by going back to the values in y.
324  vector_copy(nvar,y,yout);
325 
326  bool try_step_again=true;
327  while (try_step_again) {
328  try_step_again=false;
329 
330  if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt)) {
331  h0=dt;
332  final_step=true;
333  } else{
334  final_step=false;
335  }
336 
337  step_status=this->stepp->step(t0,h0,nvar,yout,dydx,yout,yerr,
338  dydx_out,derivs);
339 
340  // [GSL] Return if the stepper indicates a pointer or user function
341  // failure
342  if (step_status == exc_efault || step_status==exc_ebadfunc) {
343  return step_status;
344  }
345 
346  // [GSL] Check for stepper internal failure
347  if (step_status!=success) {
348 
349  // [GSL] Stepper was unable to calculate step. Try decreasing
350  // stepsize
351  double h_old=h0;
352  h0*=0.5;
353 
354  // [GSL] Check that an actual decrease in h0 occured and that
355  // the suggested h0 will change the independent variable
356  // will change by at least 1ulp
357  //double t_curr=GSL_COERCE_DBL(t);
358  //double t_next=GSL_COERCE_DBL(t+h0);
359  double t_curr=t;
360  double t_next=t+h0;
361 
362  if (fabs(h0) < fabs(h_old) && t_next != t_curr) {
363  // [GSL] Step was decreased. Undo step, and try again with
364  // new h0.
365  vector_copy(nvar,y,yout);
366  failed_steps++;
367  try_step_again=true;
368  } else {
369  // [GSL] notify user of step-size which caused the failure
370  h=h0;
371  // [GSL] restore original t value
372  t=t0;
373  return step_status;
374  }
375  }
376 
377  if (try_step_again==false) {
378 
379  count++;
380  last_step=h0;
381  if (final_step) {
382  t=t1;
383  } else{
384  t=t0+h0;
385  }
386 
387  // [GSL] Check error and attempt to adjust the step.
388  double h_old=h0;
389  const size_t hadjust_status=con.hadjust
390  (nvar,this->stepp->get_order(),yout,yerr,dydx_out,h0);
391 
392  if (hadjust_status ==
394 
395  // [GSL] Check that the reported status is correct (i.e. an actual
396  // decrease in h0 occured) and the suggested h0 will change
397  // the independent variable by at least 1 ulp
398 
399  //double t_curr=GSL_COERCE_DBL (*t);
400  //double t_next=GSL_COERCE_DBL ((*t) + h0);
401  double t_curr=t;
402  double t_next=t+h0;
403 
404  if (fabs(h0) < fabs(h_old) && t_next != t_curr) {
405  // [GSL] Step was decreased. Undo step, and try again with
406  // new h0.
407  vector_copy(nvar,y,yout);
408  failed_steps++;
409  try_step_again=true;
410  } else {
411  /* [GSL] Can not obtain required error tolerance, and can
412  not decrease step-size any further, so give up and return
413  gsl_failure.
414  */
415  // [GSL] notify user of step-size which caused the failure
416  h=h0;
417  return gsl_failure;
418  }
419  }
420  }
421 
422  }
423 
424  if (final_step==false) {
425  // [GSL] Suggest step size for next time-step. Change of step
426  // size is not suggested in the final step, because that step
427  // can be very small compared to previous step, to reach t1.
428  h=h0;
429  }
430 
431  // This used to be return step_status, but the code never
432  // reaches this point if step_status is anything other than zero.
433  return step_status;
434  }
435 
436 #endif
437 
438  public:
439 
440  astep_gsl() {
441  this->verbose=0;
442  msize=0;
443  }
444 
445  virtual ~astep_gsl() {
446  }
447 
448  /// Control specification
450 
451  /** \brief Make an adaptive integration step of the system
452  \c derivs
453 
454  This attempts to take a step of size \c h from the point \c
455  x of an \c n-dimensional system \c derivs starting with \c
456  y. On exit, \c x and \c y contain the new values at the end
457  of the step, \c h contains the size of the step, \c dydx_out
458  contains the derivative at the end of the step, and \c yerr
459  contains the estimated error at the end of the step.
460  */
461  virtual int astep(double &x, double xmax, double &h,
462  size_t n, vec_y_t &y, vec_dydx_t &dydx_out,
463  vec_yerr_t &yerr, func_t &derivs) {
464  count=0;
465  failed_steps=0;
466  last_step=0.0;
467 
468  if (n!=msize) {
469  yout_int.resize(n);
470  dydx_int.resize(n);
471  msize=n;
472  }
473 
474  int ret=derivs(x,n,y,dydx_int);
475  if (ret!=0) return ret;
476 
477  ret=evolve_apply(x,xmax,x,h,n,y,dydx_int,yout_int,yerr,dydx_out,derivs);
478 
479  if (ret==success) {
480  for(size_t i=0;i<n;i++) {
481  y[i]=yout_int[i];
482  }
483  }
484 
485  if (this->verbose>0) {
486  std::cout << "astep_gsl step:";
487  std::cout << x << " ";
488  for(size_t j=0;j<n;j++) std::cout << y[j] << " ";
489  std::cout << std::endl;
490  }
491 
492  return ret;
493  }
494 
495  /** \brief Make an adaptive integration step of the system
496  \c derivs with derivatives
497 
498  This attempts to take a step of size \c h from the point \c x
499  of an \c n-dimensional system \c derivs starting with \c y and
500  given the initial derivatives \c dydx. On exit, \c x, \c y and
501  \c dydx contain the new values at the end of the step, \c h
502  contains the size of the step, \c dydx contains the derivative
503  at the end of the step, and \c yerr contains the estimated
504  error at the end of the step.
505  */
506  virtual int astep_derivs(double &x, double xmax, double &h,
507  size_t n, vec_y_t &y, vec_dydx_t &dydx,
508  vec_yerr_t &yerr, func_t &derivs) {
509  count=0;
510  failed_steps=0;
511  last_step=0.0;
512 
513  if (n!=msize) {
514  yout_int.resize(n);
515  dydx_int.resize(n);
516  msize=n;
517  }
518 
519  int ret=evolve_apply(x,xmax,x,h,n,y,dydx,yout_int,yerr,dydx_int,derivs);
520 
521  if (ret==success) {
522  for(size_t i=0;i<n;i++) {
523  y[i]=yout_int[i];
524  dydx[i]=dydx_int[i];
525  }
526  }
527 
528  if (this->verbose>0) {
529  std::cout << "astep_gsl step:";
530  std::cout << x << " ";
531  for(size_t j=0;j<n;j++) std::cout << y[j] << " ";
532  std::cout << std::endl;
533  }
534 
535  return ret;
536  }
537 
538  /** \brief Make an adaptive integration step of the system
539  \c derivs
540 
541  This function performs an adaptive integration step with the
542  \c n-dimensional system \c derivs and parameter \c pa. It
543  Begins at \c x with initial stepsize \c h, ensuring that the
544  step goes no farther than \c xmax. At the end of the step, the
545  size of the step taken is \c h and the new value of \c x is in
546  \c x_out. Initially, the function values and derivatives
547  should be specified in \c y and \c dydx. The function values,
548  derivatives, and the error at the end of the step are given in
549  \c yout, \c yerr, and \c dydx_out. Unlike in \c ode_step
550  objects, the objects \c y, \c yout, \c dydx, and \c dydx_out
551  must all be distinct.
552 
553  This adaptive stepper function is faster than astep() or
554  astep_derivs() because it does not require any copying of
555  vectors.
556  */
557  virtual int astep_full(double x, double xmax, double &x_out, double &h,
558  size_t n, vec_y_t &y, vec_dydx_t &dydx,
559  vec_y_t &yout, vec_yerr_t &yerr,
560  vec_dydx_t &dydx_out, func_t &derivs) {
561 
562  count=0;
563  failed_steps=0;
564  last_step=0.0;
565 
566  int ret=evolve_apply(x,xmax,x_out,h,n,y,dydx,yout,yerr,dydx_out,
567  derivs);
568 
569  if (this->verbose>0) {
570  std::cout << "astep_gsl step:";
571  std::cout << x_out << " ";
572  for(size_t j=0;j<n;j++) std::cout << yout[j] << " ";
573  std::cout << std::endl;
574  }
575 
576  return ret;
577  }
578 
579  };
580 
581 #ifndef DOXYGEN_NO_O2NS
582 }
583 #endif
584 
585 #endif
invalid pointer
Definition: err_hnd.h:57
unsigned long int count
The number of steps.
Definition: astep_gsl.h:289
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
double a_dydt
Derivative scaling factor (default 0)
Definition: astep_gsl.h:160
invalid argument supplied by user
Definition: err_hnd.h:59
double last_step
The size of the last step.
Definition: astep_gsl.h:286
Control structure for astep_gsl.
Definition: astep_gsl.h:130
virtual int astep_derivs(double &x, double xmax, double &h, size_t n, vec_y_t &y, vec_dydx_t &dydx, vec_yerr_t &yerr, func_t &derivs)
Make an adaptive integration step of the system derivs with derivatives.
Definition: astep_gsl.h:506
static const size_t hadj_dec
Recommend step decrease.
Definition: astep_gsl.h:148
unsigned long int failed_steps
The number of failed steps.
Definition: astep_gsl.h:292
bool standard
Use standard or scaled algorithm (default true)
Definition: astep_gsl.h:162
double a_y
Function scaling factor (default 1)
Definition: astep_gsl.h:158
ode_control_gsl< vec_y_t, vec_dydx_t, vec_yerr_t > con
Control specification.
Definition: astep_gsl.h:449
size_t msize
The size of the allocated vectors.
Definition: astep_gsl.h:295
int evolve_apply(double t0, double t1, double &t, double &h, size_t nvar, vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out, func_t &derivs)
Apply the evolution for the next adaptive step.
Definition: astep_gsl.h:304
virtual int astep_full(double x, double xmax, double &x_out, 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)
Make an adaptive integration step of the system derivs.
Definition: astep_gsl.h:557
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
static const size_t hadj_nil
No adjustment required.
Definition: astep_gsl.h:146
ubvector scale_abs
Scalings.
Definition: astep_gsl.h:139
static const size_t hadj_inc
Recommend step increase.
Definition: astep_gsl.h:150
Failure.
Definition: err_hnd.h:49
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Adaptive stepper [abstract base].
Definition: astep.h:49
virtual int astep(double &x, double xmax, double &h, size_t n, vec_y_t &y, vec_dydx_t &dydx_out, vec_yerr_t &yerr, func_t &derivs)
Make an adaptive integration step of the system derivs.
Definition: astep_gsl.h:461
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
vec_dydx_t dydx_int
Internal storage for dydx.
Definition: astep_gsl.h:283
double eps_rel
Relative precision (default 0)
Definition: astep_gsl.h:156
vec_y_t yout_int
Temporary storage for yout.
Definition: astep_gsl.h:280
double eps_abs
Absolute precision (default )
Definition: astep_gsl.h:154
problem with user-supplied function
Definition: err_hnd.h:69
int set_scale(size_t nscal, const svec_t &scale)
Set the scaling for each differential equation.
Definition: astep_gsl.h:177
Success.
Definition: err_hnd.h:47
Adaptive ODE stepper (GSL)
Definition: astep_gsl.h:268

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