ode_it_solve.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 #ifndef O2SCL_ODE_IT_SOLVE_H
24 #define O2SCL_ODE_IT_SOLVE_H
25 
26 /** \file ode_it_solve.h
27  \brief File defining \ref o2scl::ode_it_solve
28 */
29 
30 #include <iostream>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/vector_proxy.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 
37 #include <o2scl/misc.h>
38 #include <o2scl/test_mgr.h>
39 #include <o2scl/linear_solver.h>
40 #include <o2scl/vector.h>
41 
42 #ifndef DOXYGEN_NO_O2NS
43 namespace o2scl {
44 #endif
45 
46  /// Function for iterative solving of ODEs
47  typedef std::function<double
50 
51  /** \brief ODE solver using a generic linear solver to solve
52  finite-difference equations
53 
54  \future Set up convergence error if it goes beyond max iterations
55  \future Create a GSL-like set() and iterate() interface
56  \future Implement as a child of ode_bv_solve ?
57  \future Max and average tolerance?
58  \future Allow the user to ensure that the solver doesn't
59  apply the full correction
60  */
61  template <class func_t=ode_it_funct11,
64  class matrix_row_t=boost::numeric::ublas::matrix_row
66  class solver_vec_t=boost::numeric::ublas::vector<double> ,
67  class solver_mat_t=boost::numeric::ublas::matrix<double> >
68  class ode_it_solve {
69 
70  public:
71 
72  bool make_mats;
73 
74  ode_it_solve() {
75  eps_rel=sqrt(std::numeric_limits<double>::epsilon());
76  eps_min=1.0e-15;
77  niter=30;
78  tol_rel=1.0e-8;
79  verbose=0;
81  alpha=1.0;
82  make_mats=false;
83  fd=0;
84  fl=0;
85  fr=0;
86  }
87 
88  virtual ~ode_it_solve() {}
89 
90  /// Set level of output (default 0)
91  int verbose;
92 
93  /** \brief Stepsize for finite differencing (default \f$ 10^{-4} \f$)
94  */
95  double eps_rel;
96 
97  /** \brief Minimum stepsize for finite differencing (default \f$
98  10^{-15} \f$)
99  */
100  double eps_min;
101 
102  /// Tolerance (default \f$ 10^{-8} \f$)
103  double tol_rel;
104 
105  /// Maximum number of iterations (default 30)
106  size_t niter;
107 
108  /// Set the linear solver
110  solver=&ls;
111  return 0;
112  }
113 
114  /// Size of correction to apply (default 1.0)
115  double alpha;
116 
117  /** \brief Solve \c derivs with boundary conditions \c left and
118  \c right
119 
120  Given a grid of size \c n_grid and \c n_eq differential equations,
121  solve them by relaxation. The grid is specified in \c x, which
122  is a vector of size \c n_grid. The differential equations are
123  given in \c derivs, the boundary conditions on the left hand
124  side in \c left, and the boundary conditions on the right hand
125  side in \c right. The number of boundary conditions on the left
126  hand side is \c nb_left, and the number of boundary conditions on
127  the right hand side should be <tt>n_eq-nb_left</tt>. The initial
128  guess for the solution, a matrix of size <tt>[n_grid][n_eq]</tt>
129  should be given in \c y. Upon success, \c y will contain an
130  approximate solution of the differential equations. The matrix
131  \c mat is workspace of size <tt>[n_grid*n_eq][n_grid*n_eq]</tt>, and
132  the vectors \c rhs and \c y are workspace of size
133  <tt>[n_grid*n_eq]</tt>.
134  */
135  int solve(size_t n_grid, size_t n_eq, size_t nb_left, vec_t &x,
136  mat_t &y, func_t &derivs, func_t &left, func_t &right,
137  solver_mat_t &mat, solver_vec_t &rhs, solver_vec_t &dy) {
138 
139  // Store the functions for simple derivatives
140  fd=&derivs;
141  fl=&left;
142  fr=&right;
143 
144  /// Function derivatives for iterative solving of ODEs
145  typedef std::function<double
146  (size_t,size_t,double,matrix_row_t &)> ode_it_dfunct11;
147 
148  ode_it_dfunct11 d2_derivs=std::bind
149  (std::mem_fn<double(size_t,size_t,double,matrix_row_t &)>
150  (&ode_it_solve::fd_derivs),this,std::placeholders::_1,
151  std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
152  ode_it_dfunct11 d2_left=std::bind
153  (std::mem_fn<double(size_t,size_t,double,matrix_row_t &)>
154  (&ode_it_solve::fd_left),this,std::placeholders::_1,
155  std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
156  ode_it_dfunct11 d2_right=std::bind
157  (std::mem_fn<double(size_t,size_t,double,matrix_row_t &)>
158  (&ode_it_solve::fd_right),this,std::placeholders::_1,
159  std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
160 
161  return solve_derivs(n_grid,n_eq,nb_left,x,y,derivs,left,right,
162  d2_derivs,d2_left,d2_right,mat,rhs,dy);
163  }
164 
165  /** \brief Solve \c derivs with boundary conditions \c left and
166  \c right
167 
168  Given a grid of size \c n_grid and \c n_eq differential equations,
169  solve them by relaxation. The grid is specified in \c x, which
170  is a vector of size \c n_grid. The differential equations are
171  given in \c derivs, the boundary conditions on the left hand
172  side in \c left, and the boundary conditions on the right hand
173  side in \c right. The number of boundary conditions on the left
174  hand side is \c nb_left, and the number of boundary conditions on
175  the right hand side should be <tt>n_eq-nb_left</tt>. The initial
176  guess for the solution, a matrix of size <tt>[n_grid][n_eq]</tt>
177  should be given in \c y. Upon success, \c y will contain an
178  approximate solution of the differential equations. The matrix
179  \c mat is workspace of size <tt>[n_grid*n_eq][n_grid*n_eq]</tt>, and
180  the vectors \c rhs and \c y are workspace of size
181  <tt>[n_grid*n_eq]</tt>.
182  */
183  template<class dfunc_t>
184  int solve_derivs(size_t n_grid, size_t n_eq, size_t nb_left, vec_t &x,
185  mat_t &y, func_t &derivs, func_t &left, func_t &right,
186  dfunc_t &d_derivs, dfunc_t &d_left, dfunc_t &d_right,
187  solver_mat_t &mat, solver_vec_t &rhs, solver_vec_t &dy) {
188 
189  // Variable index
190  size_t ix;
191 
192  // Number of RHS boundary conditions
193  size_t nb_right=n_eq-nb_left;
194 
195  // Number of variables
196  size_t nvars=n_grid*n_eq;
197 
198  bool done=false;
199  for(size_t it=0;done==false && it<niter;it++) {
200 
201  ix=0;
202 
203  for(size_t i=0;i<nvars;i++) {
204  for(size_t j=0;j<nvars;j++) {
205  mat(i,j)=0.0;
206  }
207  }
208 
209  // Construct the entries corresponding to the LHS boundary.
210  // This makes the first nb_left rows of the matrix.
211  for(size_t i=0;i<nb_left;i++) {
212  matrix_row_t yk=o2scl::matrix_row<mat_t,matrix_row_t>(y,0);
213  rhs[ix]=-left(i,x[0],yk);
214  for(size_t j=0;j<n_eq;j++) {
215  mat(ix,j)=d_left(i,j,x[0],yk);
216  }
217  ix++;
218  }
219 
220  // Construct the matrix entries for the internal points
221  // This loop adds n_grid-1 sets of n_eq rows
222  for(size_t k=0;k<n_grid-1;k++) {
223  size_t kp1=k+1;
224  double tx=(x[kp1]+x[k])/2.0;
225  double dx=x[kp1]-x[k];
226  matrix_row_t yk=o2scl::matrix_row<mat_t,matrix_row_t>(y,k);
227  matrix_row_t ykp1=o2scl::matrix_row<mat_t,matrix_row_t>(y,k+1);
228 
229  for(size_t i=0;i<n_eq;i++) {
230 
231  rhs[ix]=y(k,i)-y(kp1,i)+(x[kp1]-x[k])*
232  (derivs(i,tx,ykp1)+derivs(i,tx,yk))/2.0;
233 
234  size_t lhs=k*n_eq;
235  for(size_t j=0;j<n_eq;j++) {
236  mat(ix,lhs+j)=-d_derivs(i,j,tx,yk)*dx/2.0;
237  mat(ix,lhs+j+n_eq)=-d_derivs(i,j,tx,ykp1)*dx/2.0;
238  if (i==j) {
239  mat(ix,lhs+j)=mat(ix,lhs+j)-1.0;
240  mat(ix,lhs+j+n_eq)=mat(ix,lhs+j+n_eq)+1.0;
241  }
242  }
243 
244  ix++;
245 
246  }
247 
248  }
249 
250  // Construct the entries corresponding to the RHS boundary
251  // This makes the last nb_right rows of the matrix.
252  for(size_t i=0;i<nb_right;i++) {
253  matrix_row_t ylast=o2scl::matrix_row<mat_t,matrix_row_t>(y,n_grid-1);
254  size_t lhs=n_eq*(n_grid-1);
255 
256  rhs[ix]=-right(i,x[n_grid-1],ylast);
257 
258  for(size_t j=0;j<n_eq;j++) {
259  mat(ix,lhs+j)=d_right(i,j,x[n_grid-1],ylast);
260  }
261 
262  ix++;
263 
264  }
265 
266  // Compute correction by calling the linear solver
267 
268  if (verbose>3) {
269  std::cout << "Matrix: " << std::endl;
270  for(size_t i=0;i<nvars;i++) {
271  for(size_t j=0;j<nvars;j++) {
272  std::cout << mat(i,j) << " ";
273  }
274  std::cout << std::endl;
275  }
276  std::cout << "Deviations:" << std::endl;
277  for(size_t i=0;i<nvars;i++) {
278  std::cout << rhs[i] << std::endl;
279  }
280  }
281 
282  if (make_mats) return 0;
283 
284  solver->solve(ix,mat,rhs,dy);
285 
286  if (verbose>3) {
287  std::cout << "Corrections:" << std::endl;
288  for(size_t i=0;i<nvars;i++) {
289  std::cout << dy[i] << std::endl;
290  }
291  std::cout << "Press a key and press enter to continue: "
292  << std::endl;
293  char ch;
294  std::cin >> ch;
295  }
296 
297  // Apply correction and compute its size
298 
299  double res=0.0;
300  ix=0;
301 
302  for(size_t igrid=0;igrid<n_grid;igrid++) {
303  for(size_t ieq=0;ieq<n_eq;ieq++) {
304  y(igrid,ieq)+=alpha*dy[ix];
305  res+=dy[ix]*dy[ix];
306  ix++;
307  }
308  }
309 
310  if (verbose>0) {
311  // Since we're in the o2scl namespace, we explicitly
312  // specify std::sqrt() here
313  std::cout << "ode_it_solve: " << it << " " << std::sqrt(res) << " "
314  << tol_rel << std::endl;
315  if (verbose>1) {
316  char ch;
317  std::cout << "Press a key and type enter to continue. ";
318  std::cin >> ch;
319  }
320  }
321 
322  // If the correction has become small enough, we're done
323  if (std::sqrt(res)<=tol_rel) done=true;
324  }
325 
326  if (done==false) {
327  O2SCL_ERR("Exceeded number of iterations in solve().",
329  }
330 
331  return 0;
332  }
333 
334  /// Default linear solver
336 
337  protected:
338 
339  /// \name Storage for functions
340  //@{
341  func_t *fl, *fr, *fd;
342  //@}
343 
344  /// Solver
346 
347  /** \brief Compute the derivatives of the LHS boundary conditions
348 
349  This function computes \f$ \partial f_{left,\mathrm{ieq}} / \partial
350  y_{\mathrm{ivar}} \f$
351  */
352  virtual double fd_left(size_t ieq, size_t ivar, double x, matrix_row_t &y) {
353 
354  double ret, dydx;
355 
356  double h=eps_rel*fabs(y[ivar]);
357  if (h<eps_min) h=eps_min;
358  if (h==0.0) h=eps_rel;
359 
360  y[ivar]+=h;
361  ret=(*fl)(ieq,x,y);
362 
363  y[ivar]-=h;
364  ret-=(*fl)(ieq,x,y);
365 
366  ret/=h;
367  return ret;
368  }
369 
370  /** \brief Compute the derivatives of the RHS boundary conditions
371 
372  This function computes \f$ \partial f_{right,\mathrm{ieq}} / \partial
373  y_{\mathrm{ivar}} \f$
374  */
375  virtual double fd_right(size_t ieq, size_t ivar, double x, matrix_row_t &y) {
376 
377  double ret, dydx;
378 
379  double h=eps_rel*fabs(y[ivar]);
380  if (h<eps_min) h=eps_min;
381  if (h==0.0) h=eps_rel;
382 
383  y[ivar]+=h;
384  ret=(*fr)(ieq,x,y);
385 
386  y[ivar]-=h;
387  ret-=(*fr)(ieq,x,y);
388 
389  ret/=h;
390  return ret;
391  }
392 
393  /** \brief Compute the finite-differenced part of the
394  differential equations
395 
396  This function computes \f$ \partial f_{\mathrm{ieq}} / \partial
397  y_{\mathrm{ivar}} \f$
398  */
399  virtual double fd_derivs(size_t ieq, size_t ivar, double x,
400  matrix_row_t &y) {
401 
402  double ret, dydx;
403 
404  double h=eps_rel*fabs(y[ivar]);
405  if (h<eps_min) h=eps_min;
406  if (h==0.0) h=eps_rel;
407 
408  y[ivar]+=h;
409  ret=(*fd)(ieq,x,y);
410 
411  y[ivar]-=h;
412  ret-=(*fd)(ieq,x,y);
413 
414  ret/=h;
415 
416  return ret;
417  }
418 
419  };
420 
421 #ifndef DOXYGEN_NO_O2NS
422 }
423 #endif
424 
425 #endif
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
ODE solver using a generic linear solver to solve finite-difference equations.
Definition: ode_it_solve.h:68
int set_solver(o2scl_linalg::linear_solver< solver_vec_t, solver_mat_t > &ls)
Set the linear solver.
Definition: ode_it_solve.h:109
double alpha
Size of correction to apply (default 1.0)
Definition: ode_it_solve.h:115
exceeded max number of iterations
Definition: err_hnd.h:73
int verbose
Set level of output (default 0)
Definition: ode_it_solve.h:91
o2scl_linalg::linear_solver< solver_vec_t, solver_mat_t > * solver
Solver.
Definition: ode_it_solve.h:345
virtual double fd_left(size_t ieq, size_t ivar, double x, matrix_row_t &y)
Compute the derivatives of the LHS boundary conditions.
Definition: ode_it_solve.h:352
double eps_rel
Stepsize for finite differencing (default )
Definition: ode_it_solve.h:95
int solve(size_t n_grid, size_t n_eq, size_t nb_left, vec_t &x, mat_t &y, func_t &derivs, func_t &left, func_t &right, solver_mat_t &mat, solver_vec_t &rhs, solver_vec_t &dy)
Solve derivs with boundary conditions left and right.
Definition: ode_it_solve.h:135
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int solve_derivs(size_t n_grid, size_t n_eq, size_t nb_left, vec_t &x, mat_t &y, func_t &derivs, func_t &left, func_t &right, dfunc_t &d_derivs, dfunc_t &d_left, dfunc_t &d_right, solver_mat_t &mat, solver_vec_t &rhs, solver_vec_t &dy)
Solve derivs with boundary conditions left and right.
Definition: ode_it_solve.h:184
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
Definition: vector.h:2221
virtual double fd_derivs(size_t ieq, size_t ivar, double x, matrix_row_t &y)
Compute the finite-differenced part of the differential equations.
Definition: ode_it_solve.h:399
virtual void solve(size_t n, mat_t &a, vec_t &b, vec_t &x)=0
Solve square linear system of size n.
std::function< double(size_t, double, boost::numeric::ublas::matrix_row< boost::numeric::ublas::matrix< double > > &)> ode_it_funct11
Function for iterative solving of ODEs.
Definition: ode_it_solve.h:49
double tol_rel
Tolerance (default )
Definition: ode_it_solve.h:103
double eps_min
Minimum stepsize for finite differencing (default )
Definition: ode_it_solve.h:100
virtual double fd_right(size_t ieq, size_t ivar, double x, matrix_row_t &y)
Compute the derivatives of the RHS boundary conditions.
Definition: ode_it_solve.h:375
o2scl_linalg::linear_solver_HH< solver_vec_t, solver_mat_t > def_solver
Default linear solver.
Definition: ode_it_solve.h:335
size_t niter
Maximum number of iterations (default 30)
Definition: ode_it_solve.h:106

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