ode_bsimp_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/bsimp.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 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,
40  * Boston, MA 02110-1301, USA.
41  */
42 
43 #ifndef O2SCL_GSL_BSIMP_H
44 #define O2SCL_GSL_BSIMP_H
45 
46 /** \file ode_bsimp_gsl.h
47  \brief File defining \ref o2scl::ode_bsimp_gsl
48 */
49 
50 #include <gsl/gsl_math.h>
51 
52 #include <o2scl/err_hnd.h>
53 #include <o2scl/ode_step.h>
54 #include <o2scl/ode_jac_funct.h>
55 #include <o2scl/lu.h>
56 #include <o2scl/permutation.h>
57 
58 #ifndef DOXYGEN_NO_O2NS
59 namespace o2scl {
60 #endif
61 
62  /** \brief Bulirsch-Stoer implicit ODE stepper (GSL)
63 
64  Implicit Bulirsch-Stoer method of Bader and Deuflhard (\ref
65  Bader83).
66 
67  \note The variable <tt>h_next</tt> was defined in the original
68  GSL version has been removed here, as it was unused by the
69  stepper routine.
70 
71  \note At the moment, this class retains the GSL approach to
72  handling non-integer return values in the user-specified
73  derivative function. If the user-specified function returns \c
74  exc_efailed, then the stepper attempts to decrease the stepsize
75  to continue. If the user-specified function returns a non-zero
76  value other than \c exc_efailed, or if the Jacobian evaluation
77  returns any non-zero value, then the stepper aborts and returns
78  the error value without calling the error handler. This behavior
79  may change in the future.
80 
81  There is an example for the usage of this class in
82  <tt>examples/ex_stiff.cpp</tt> documented in the \ref
83  ex_stiff_sect section.
84 
85  \future More detailed documentation about the algorithm
86 
87  \future Figure out if this should be a child of ode_step or
88  astep. The function step_local() is actually its own ODE
89  stepper and could be reimplemented as an object of type ode_step.
90 
91  \future I don't like setting yerr to GSL_POSINF, there should
92  be a better way to force an adaptive stepper which is calling
93  this stepper to readjust the stepsize.
94 
95  \future The functions deuf_kchoice() and compute_weights() can
96  be moved out of this header file.
97 
98  \future Rework internal arrays
99 
100  \future Rework linear solver to be amenable to using
101  a sparse matrix solver
102  */
103  template<class func_t=ode_funct, class jac_func_t=ode_jac_funct,
106 
107  public:
108 
109  typedef boost::numeric::ublas::unbounded_array<double> ubarr;
112 
113 #ifndef DOXYGEN_INTERAL
114 
115  protected:
116 
117  /// Size of allocated vectors
118  size_t dim;
119 
120  /// Function specifying derivatives
121  func_t *funcp;
122 
123  /// Jacobian
124  jac_func_t *jfuncp;
125 
126  /// Number of entries in the Bader-Deuflhard extrapolation sequence
127  static const int sequence_count=8;
128 
129  /// Largest index of the Bader-Deuflhard extrapolation sequence
130  static const int sequence_max=7;
131 
132  /// Workspace for extrapolation
134 
135  /// Workspace for linear system matrix
137 
138  /** \brief Workspace for extrapolation
139 
140  (This state variable was named 'x' in GSL.)
141  */
142  double ex_wk[sequence_max];
143 
144  /// \name State info
145  //@{
146  size_t k_current;
147  size_t k_choice;
148  double eps;
149  //@}
150 
151  /// \name Workspace for extrapolation step
152  //@{
153  ubarr yp;
154  ubarr y_save;
155  ubarr yerr_save;
156  ubarr y_extrap_save;
157  vec_t y_extrap_sequence;
158  ubarr extrap_work;
159  vec_t dfdt;
160  vec_t y_temp;
161  vec_t delta_temp;
162  ubarr weight;
163  mat_t dfdy;
164  //@}
165 
166  /// \name Workspace for the basic stepper
167  //@{
168  vec_t rhs_temp;
169  ubarr delta;
170  //@}
171 
172  /// Order of last step
173  size_t order;
174 
175  /// Compute weights
176  int compute_weights(const ubarr &y, ubarr &w, size_t n) {
177  size_t i;
178  // w[i] is 1 if y[i] is zero and the absolute value of y[i]
179  // otherwise
180  for (i=0;i<n;i++) {
181  double u=fabs(y[i]);
182  w[i]=(u>0.0) ? u : 1.0;
183  }
184  return 0;
185  }
186 
187  /** \brief Calculate a choice for the "order" of the method, using the
188  Deuflhard criteria.
189 
190  Used in the allocate() function.
191  */
192  size_t deuf_kchoice(double eps2, size_t dimension) {
193 
194  const double safety_f=0.25;
195  const double small_eps=safety_f*eps2;
196 
197  double a_work[sequence_count];
198  double alpha[sequence_max][sequence_max];
199 
200  /* Bader-Deuflhard extrapolation sequence */
201  static const int bd_sequence[sequence_count]=
202  {2,6,10,14,22,34,50,70};
203 
204  int i, k;
205 
206  a_work[0]=bd_sequence[0]+1.0;
207 
208  for (k=0;k<sequence_max;k++) {
209  a_work[k+1]=a_work[k]+bd_sequence[k+1];
210  }
211 
212  for (i=0;i<sequence_max;i++) {
213  alpha[i][i]=1.0;
214  for (k=0;k<i;k++) {
215  const double tmp1=a_work[k+1]-a_work[i+1];
216  const double tmp2=(a_work[i+1]-a_work[0]+1.0)*(2*k+1);
217  alpha[k][i]=pow(small_eps,tmp1/tmp2);
218  }
219  }
220 
221  a_work[0]+=dimension;
222 
223  for (k=0;k<sequence_max;k++) {
224  a_work[k+1]=a_work[k]+bd_sequence[k+1];
225  }
226 
227  for (k=0;k<sequence_max-1;k++) {
228  if (a_work[k+2]>a_work[k+1]*alpha[k][k+1]) {
229  break;
230  }
231  }
232 
233  return k;
234  }
235 
236  /** \brief Polynomial extrapolation
237 
238  Compute the step of index <tt>i_step</tt> using polynomial
239  extrapolation to evaulate functions by fitting a polynomial
240  to estimates <tt>(x_i,y_i)</tt> and output the result to
241  <tt>y_0</tt> and <tt>y_0_err</tt>.
242 
243  The index <tt>i_step</tt> begins with zero.
244 
245  */
246  int poly_extrap(ubmatrix &dloc, const double x[],
247  const unsigned int i_step, const double x_i,
248  const vec_t &y_i, vec_t &y_0, vec_t &y_0_err,
249  ubarr &work) {
250  size_t j, k;
251 
252  o2scl::vector_copy(dim,y_i,y_0_err);
253  o2scl::vector_copy(dim,y_i,y_0);
254 
255  if (i_step == 0) {
256 
257  for (j=0;j<dim;j++) {
258  dloc(0,j)=y_i[j];
259  }
260 
261  } else {
262 
263  o2scl::vector_copy(dim,y_i,work);
264 
265  for (k=0;k<i_step;k++) {
266  double deltaloc=1.0/(x[i_step-k-1]-x_i);
267  const double f1=deltaloc*x_i;
268  const double f2=deltaloc*x[i_step-k-1];
269 
270  for (j=0;j<dim;j++) {
271  const double q_kj=dloc(k,j);
272  dloc(k,j)=y_0_err[j];
273  deltaloc=work[j]-q_kj;
274  y_0_err[j]=f1*deltaloc;
275  work[j]=f2*deltaloc;
276  y_0[j]+=y_0_err[j];
277  }
278  }
279 
280  for (j=0;j<dim;j++) {
281  dloc(i_step,j)=y_0_err[j];
282  }
283  }
284  return 0;
285  }
286 
287  /** \brief Basic implicit Bulirsch-Stoer step
288 
289  Divide the step <tt>h_total</tt> into <tt>n_step</tt> smaller
290  steps and do the Bader-Deuflhard semi-implicit iteration. This
291  function starts at <tt>t0</tt> with function values
292  <tt>y</tt>, derivatives <tt>yp_loc</tt>, and information from
293  the Jacobian to compute the final value <tt>y_out</tt>.
294  */
295  int step_local(const double t0, const double h_total,
296  const unsigned int n_step, const ubarr &y,
297  const ubarr &yp_loc, const vec_t &dfdt_loc,
298  const mat_t &dfdy_loc, vec_t &y_out) {
299 
300  const double h=h_total/n_step;
301  double t=t0+h;
302 
303  double sum;
304 
305  /* This is the factor sigma referred to in equation 3.4 of the
306  paper. A relative change in y exceeding sigma indicates a
307  runaway behavior. According to the authors suitable values
308  for sigma are >>1. I have chosen a value of 100*dim. BJG
309  */
310  const double max_sum=100.0*dim;
311 
312  int signum, status;
313  size_t i, j;
314  size_t n_inter;
315 
316  /* Calculate the matrix for the linear system. */
317  for (i=0;i<dim;i++) {
318  for (j=0;j<dim;j++) {
319  a_mat(i,j)=-h*dfdy_loc(i,j);
320  }
321  a_mat(i,i)=a_mat(i,i)+1.0;
322  }
323 
324  /* LU decomposition for the linear system. */
325 
326  o2scl::permutation p_vec(dim);
327 
328  o2scl_linalg::LU_decomp(dim,a_mat,p_vec,signum);
329 
330  /* Compute weighting factors */
331  compute_weights(y,weight,dim);
332 
333  /* Initial step. */
334 
335  for (i=0;i<dim;i++) {
336  y_temp[i]=h*(yp_loc[i]+h*dfdt_loc[i]);
337  }
338 
339  o2scl_linalg::LU_solve(dim,a_mat,p_vec,y_temp,delta_temp);
340 
341  sum=0.0;
342  for (i=0;i<dim;i++) {
343  const double di=delta_temp[i];
344  delta[i]=di;
345  y_temp[i]=y[i]+di;
346  sum+=fabs(di)/weight[i];
347  }
348  if (sum>max_sum) {
349  return exc_efailed;
350  }
351 
352  /* Intermediate steps. */
353 
354  status=(*funcp)(t,dim,y_temp,y_out);
355  if (status) {
356  return status;
357  }
358 
359  for (n_inter=1;n_inter<n_step;n_inter++) {
360 
361  for (i=0;i<dim;i++) {
362  rhs_temp[i]=h*y_out[i]-delta[i];
363  }
364 
365  o2scl_linalg::LU_solve(dim,a_mat,p_vec,rhs_temp,delta_temp);
366 
367  sum=0.0;
368  for (i=0;i<dim;i++) {
369  delta[i]+=2.0*delta_temp[i];
370  y_temp[i]+=delta[i];
371  sum+=fabs(delta[i])/weight[i];
372  }
373  if (sum>max_sum) {
374  return exc_efailed;
375  }
376 
377  t+=h;
378 
379  status=(*funcp)(t,dim,y_temp,y_out);
380  if (status) {
381  return status;
382  }
383  }
384 
385 
386  /* Final step. */
387 
388  for (i=0;i<dim;i++) {
389  rhs_temp[i]=h*y_out[i]-delta[i];
390  }
391 
392  o2scl_linalg::LU_solve(dim,a_mat,p_vec,rhs_temp,delta_temp);
393 
394  sum=0.0;
395  for (i=0;i<dim;i++) {
396  y_out[i]=y_temp[i]+delta_temp[i];
397  sum+=fabs(delta_temp[i])/weight[i];
398  }
399 
400  if (sum>max_sum) {
401  return exc_efailed;
402  }
403 
404  return success;
405  }
406 
407  /// Allocate memory for a system of size \c n
408  int allocate(size_t n) {
409 
410  if (dim>0) free();
411 
412  dim=n;
413 
414  d.resize(sequence_max,n);
415  a_mat.resize(n,n);
416 
417  yp.resize(n);
418 
419  // AWS, 12/2/08 - This was added to ensure memory reallocation
420  // resets the stepper just like reset() does
421  for(size_t i=0;i<n;i++) yp[i]=0.0;
422 
423  y_save.resize(n);
424  yerr_save.resize(n);
425  y_extrap_save.resize(n);
426  extrap_work.resize(n);
427  weight.resize(n);
428 
429  dfdy.resize(n,n);
430  dfdt.resize(n);
431  y_extrap_sequence.resize(n);
432  y_temp.resize(n);
433  rhs_temp.resize(n);
434  delta_temp.resize(n);
435 
436  delta.resize(n);
437 
438  double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
439 
440  // This choice of epsilon is not necessarily optimal, it has
441  // a "FIXME" comment in the original GSL code
442  size_t k_choice_loc=deuf_kchoice(sqrt_dbl_eps,n);
443  k_choice=k_choice_loc;
444  k_current=k_choice_loc;
445  order=2*k_choice_loc;
446 
447  return 0;
448  }
449 
450  /// Free allocated memory
451  void free() {
452  if (dim>0) {
453  delta.resize(0);
454  weight.resize(0);
455  extrap_work.resize(0);
456  y_extrap_save.resize(0);
457  y_save.resize(0);
458  yerr_save.resize(0);
459  yp.resize(0);
460  d.resize(0,0);
461  dim=0;
462  }
463  }
464 
465 #endif
466 
467  public:
468 
469  ode_bsimp_gsl() {
470  dim=0;
471  verbose=0;
472  }
473 
474  virtual ~ode_bsimp_gsl() {
475  if (dim>0) free();
476  }
477 
478  /// Verbose parameter
479  int verbose;
480 
481  /// Reset stepper
482  int reset() {
483  for(size_t i=0;i<dim;i++) yp[i]=0.0;
484  return success;
485  }
486 
487  /** \brief Perform an integration step
488 
489  Given initial value of the n-dimensional function in \c y and
490  the derivative in \c dydx (which must generally be computed
491  beforehand) at the point \c x, take a step of size \c h giving
492  the result in \c yout, the uncertainty in \c yerr, and the new
493  derivative in \c dydx_out (at \c x+h) using function \c derivs
494  to calculate derivatives. Implementations which do not
495  calculate \c yerr and/or \c dydx_out do not reference these
496  variables so that a blank \c vec_t can be given. All of the
497  implementations allow \c yout=y and \c dydx_out=dydx if
498  necessary.
499  */
500  virtual int step(double x, double h, size_t n, vec_t &y, vec_t &dydx,
501  vec_t &yout, vec_t &yerr, vec_t &dydx_out,
502  func_t &derivs, jac_func_t &jac) {
503 
504  int ret;
505 
506  funcp=&derivs;
507  jfuncp=&jac;
508 
509  if (n!=dim) allocate(n);
510 
511  /* Bader-Deuflhard extrapolation sequence */
512  static const int bd_sequence[sequence_count]={2,6,10,14,22,34,50,70};
513 
514  double t_local=x;
515 
516  size_t i, k;
517 
518  if (h+t_local==t_local) {
519  // This section is labeled with a "FIXME" comment in the
520  // original GSL code. I'm not sure why, but an error is
521  // sensible here.
522  O2SCL_ERR("Stepsize underflow in ode_bsimp_gsl::step().",
523  exc_eundrflw);
524  }
525 
526  /* Save inputs */
527  o2scl::vector_copy(dim,y,y_extrap_save);
528  o2scl::vector_copy(dim,y,y_save);
529  o2scl::vector_copy(dim,yerr,yerr_save);
530 
531  // Copy derivative
532  o2scl::vector_copy(dim,dydx,yp);
533 
534  // Evaluate the Jacobian for the system. */
535  ret=jac(t_local,dim,y,dfdy,dfdt);
536  if (ret!=success) {
537  return ret;
538  }
539 
540  /* Make a series of refined extrapolations, up to the specified
541  maximum order, which was calculated based on the Deuflhard
542  criterion in the deuf_kchoice() function (which is called by
543  allocate() ).
544  */
545  for (k=0;k<=k_current;k++) {
546 
547  const unsigned int N=bd_sequence[k];
548  const double r=(h/N);
549  const double x_k=r*r;
550 
551  // Each step computes a value of y_extrap_sequence,
552  // using the number of sub-steps dictated by
553  // the BD sequence
554  int status=step_local(t_local,h,N,y_extrap_save,yp,
555  dfdt,dfdy,y_extrap_sequence);
556 
557  if (verbose>=2) {
558  std::cout << "ode_bsimp_gsl: " << k << "/" << k_current << " "
559  << "status=" << status << std::endl;
560  }
561 
562  if (status==exc_efailed) {
563  /* If the local step fails, set the error to infinity in
564  order to force a reduction in the step size
565  */
566  for (i=0;i<dim;i++) {
567  yerr[i]=GSL_POSINF;
568  }
569  break;
570  } else if (status!=success) {
571  return status;
572  }
573 
574  // Use the information in y_extrap_sequence to compute
575  // the new value of y and yerr .
576  ex_wk[k]=x_k;
577  poly_extrap(d,ex_wk,k,x_k,y_extrap_sequence,y,yerr,extrap_work);
578  }
579 
580  /* Evaluate dydt_out[]. */
581 
582  ret=derivs(t_local+h,dim,y,dydx_out);
583 
584  // If we failed, copy the old values back to y and yerr
585  if (ret!=success) {
586  o2scl::vector_copy(dim,y_save,y);
587  o2scl::vector_copy(dim,yerr_save,yerr);
588  return ret;
589  }
590 
591  return success;
592  }
593 
594  };
595 
596 #ifndef DOXYGEN_NO_O2NS
597 }
598 #endif
599 
600 #endif
boost::numeric::ublas::matrix< double >
o2scl::ode_bsimp_gsl::funcp
func_t * funcp
Function specifying derivatives.
Definition: ode_bsimp_gsl.h:121
o2scl::ode_bsimp_gsl
Bulirsch-Stoer implicit ODE stepper (GSL)
Definition: ode_bsimp_gsl.h:105
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl::ode_bsimp_gsl::sequence_max
static const int sequence_max
Largest index of the Bader-Deuflhard extrapolation sequence.
Definition: ode_bsimp_gsl.h:130
boost::numeric::ublas::vector< double >
o2scl::ode_bsimp_gsl::step
virtual int step(double x, double h, size_t n, vec_t &y, vec_t &dydx, vec_t &yout, vec_t &yerr, vec_t &dydx_out, func_t &derivs, jac_func_t &jac)
Perform an integration step.
Definition: ode_bsimp_gsl.h:500
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
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_bsimp_gsl::free
void free()
Free allocated memory.
Definition: ode_bsimp_gsl.h:451
o2scl_linalg::LU_decomp
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
o2scl::ode_bsimp_gsl::a_mat
ubmatrix a_mat
Workspace for linear system matrix.
Definition: ode_bsimp_gsl.h:136
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::ode_bsimp_gsl::sequence_count
static const int sequence_count
Number of entries in the Bader-Deuflhard extrapolation sequence.
Definition: ode_bsimp_gsl.h:127
o2scl::ode_bsimp_gsl::compute_weights
int compute_weights(const ubarr &y, ubarr &w, size_t n)
Compute weights.
Definition: ode_bsimp_gsl.h:176
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::exc_eundrflw
@ exc_eundrflw
underflow
Definition: err_hnd.h:81
o2scl::ode_bsimp_gsl::reset
int reset()
Reset stepper.
Definition: ode_bsimp_gsl.h:482
o2scl::ode_funct
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
Ordinary differential equation function in src/ode/ode_funct.h.
Definition: ode_funct.h:46
o2scl::ode_bsimp_gsl::dim
size_t dim
Size of allocated vectors.
Definition: ode_bsimp_gsl.h:118
o2scl::ode_bsimp_gsl::poly_extrap
int poly_extrap(ubmatrix &dloc, const double x[], const unsigned int i_step, const double x_i, const vec_t &y_i, vec_t &y_0, vec_t &y_0_err, ubarr &work)
Polynomial extrapolation.
Definition: ode_bsimp_gsl.h:246
o2scl::ode_bsimp_gsl::d
ubmatrix d
Workspace for extrapolation.
Definition: ode_bsimp_gsl.h:133
o2scl::ode_bsimp_gsl::order
size_t order
Order of last step.
Definition: ode_bsimp_gsl.h:173
o2scl::ode_bsimp_gsl::step_local
int step_local(const double t0, const double h_total, const unsigned int n_step, const ubarr &y, const ubarr &yp_loc, const vec_t &dfdt_loc, const mat_t &dfdy_loc, vec_t &y_out)
Basic implicit Bulirsch-Stoer step.
Definition: ode_bsimp_gsl.h:295
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::ode_bsimp_gsl::ex_wk
double ex_wk[sequence_max]
Workspace for extrapolation.
Definition: ode_bsimp_gsl.h:142
o2scl::ode_jac_funct
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &, boost::numeric::ublas::vector< double > &) > ode_jac_funct
Functor for ODE Jacobians in src/ode/ode_jac_funct.h.
Definition: ode_jac_funct.h:40
o2scl::ode_bsimp_gsl::verbose
int verbose
Verbose parameter.
Definition: ode_bsimp_gsl.h:479
o2scl_linalg::LU_solve
int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
Solve a linear system after LU decomposition.
Definition: lu_base.h:335
o2scl::ode_bsimp_gsl::deuf_kchoice
size_t deuf_kchoice(double eps2, size_t dimension)
Calculate a choice for the "order" of the method, using the Deuflhard criteria.
Definition: ode_bsimp_gsl.h:192
o2scl::ode_bsimp_gsl::jfuncp
jac_func_t * jfuncp
Jacobian.
Definition: ode_bsimp_gsl.h:124
o2scl::ode_bsimp_gsl::allocate
int allocate(size_t n)
Allocate memory for a system of size n.
Definition: ode_bsimp_gsl.h:408

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