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

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