ode_bv_multishoot.h
Go to the documentation of this file.
1  /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2008-2020, Julien Garaud and 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_BV_MULTISHOOT_H
24 #define O2SCL_ODE_BV_MULTISHOOT_H
25 
26 /** \file ode_bv_multishoot.h
27  \brief File defining \ref o2scl::ode_bv_multishoot
28 */
29 
30 #include <string>
31 #include <o2scl/astep.h>
32 #include <o2scl/astep_gsl.h>
33 #include <o2scl/ode_iv_solve.h>
34 #include <o2scl/gsl_mroot_hybrids.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Solve a ODE boundary value problem by multishooting
41 
42  This class is experimental.
43 
44  \todo Improve documentation a little and create
45  testing code
46  */
47  template <class func_t=ode_funct<>,
48  class vec_t=ubvector, class alloc_vec_t=ubvector,
49  class alloc_t=ubvector_alloc,class vec_int_t=ubvector_int_base,
50  class mat_t=ubmatrix> class ode_bv_multishoot {
51 
52  public:
53 
55  oisp=&def_ois;
56  mrootp=&def_mroot;
57  }
58 
59  virtual ~ode_bv_multishoot() {
60  }
61 
62  /* Solve the boundary value problem on a mesh */
63  virtual int solve(vec_t &mesh, int &n_func, vec_t &y_start,
64  func_t &left_b, func_t &right_b, func_t &extra_b,
65  func_t &derivs, vec_t &x_save, mat_t &y_save) {
66 
67  /* Make copies of the input for later access */
68  this->l_mesh=&mesh;
69  this->l_n_func=&n_func;
70  this->l_y_start=&y_start;
71  this->l_left_b=&left_b;
72  this->l_right_b=&right_b;
73  this->l_extra_b=&extra_b;
74  this->l_derivs=&derivs;
75  this->l_x_save=&x_save;
76  this->l_y_save=&y_save;
77 
78  /* vector of variables */
79  int nsolve=y_start.size();
80  ubvector sx(nsolve),sy(nsolve);
81  sx=y_start;
82 
83  /* Equation solver */
84  mm_funct_mfptr<ode_bv_multishoot<func_t,vec_t,
85  alloc_vec_t,alloc_t,vec_int_t> >
86  mfm(this,&ode_bv_multishoot<func_t,vec_t,
87  alloc_vec_t,alloc_t,vec_int_t>::solve_fun);
88 
89  /* Run multishooting and save at the last step */
90  this->save=false;
91  int ret=this->mrootp->msolve(nsolve,sx,mfm);
92  this->save=true;
93  solve_fun(nsolve,sx,sy);
94 
95  return ret;
96  }
97 
98  /* Set initial value solver */
100  oisp=&ois;
101  return 0;
102  }
103 
104  /* Set equation solver */
105  int set_mroot(mroot<mm_funct<> > &root) {
106  mrootp=&root;
107  return 0;
108  }
109 
110  /* Default initial value solver */
112 
113  /* Default equation solver */
114  gsl_mroot_hybrids<mm_funct<> > def_mroot;
115 
116 #ifndef DOXYGEN_INTERNAL
117 
118  protected:
119 
120  /// The initial value solver
122 
123  /// The equation solver
124  gsl_mroot_hybrids<mm_funct<> > *mrootp;
125 
126  /// Desc
127  vec_t *l_mesh;
128  /// Desc
129  vec_t *l_y_start;
130  /// Desc
131  func_t *l_left_b;
132  /// Desc
133  func_t *l_right_b;
134  /// Desc
135  func_t *l_extra_b;
136  /// Desc
137  func_t *l_derivs;
138  /// Desc
139  int *l_n_func;
140  /// Desc
141  vec_t *l_x_save;
142  /// Desc
143  mat_t *l_y_save;
144  /// Desc
145  bool save;
146 
147  /// Function to solve
148  int solve_fun(size_t nv, const vec_t &sx, vec_t &sy) {
149 
150  double xa,xb=0.0,h;
151  ubvector y((*this->l_n_func)),y2((*this->l_n_func));
152 
153  /* We update y_start in order that derivs knows
154  all the values of parameters */
155  for(size_t i=0;i<(*this->l_y_start).size();i++) {
156  (*this->l_y_start)[i]=sx[i];
157  }
158 
159  /* A loop on each subinterval */
160  for(size_t k=0;k<(*this->l_mesh).size()-1;k++) {
161 
162  xa=(*this->l_mesh)[k];
163  xb=(*this->l_mesh)[k+1];
164  h=(xb-xa)/100.0;
165 
166  /* We load function's value at the left point of the sub-interval */
167  if (k==0) {
168  (*this->l_left_b)(xa,(*this->l_n_func),sx,y);
169  } else {
170  for(int i=0;i<(*this->l_n_func);i++)
171  y[i]=sx[i+(*this->l_n_func)*(1+k)];
172  }
173 
174  if (this->save) {
175 
176  /* iv_solver if we save */
177  int ngrid=((*this->l_x_save).size()-1)/((*this->l_mesh).size()-1)+1;
178  ubvector xxsave(ngrid);
179  ubmatrix yysave(ngrid,(*this->l_n_func));
180 
181  if (k!=((*this->l_mesh).size()-2)) {
182  xb=(*this->l_mesh)[k+1]-((*this->l_mesh)[k+1]-
183  (*this->l_mesh)[k])/ngrid;
184  }
185 
186  this->oisp->solve_grid(xa,xb,h,(*this->l_n_func),y,ngrid,
187  xxsave,yysave,(*this->l_derivs));
188 
189  for(int i=0;i<ngrid;i++) {
190  (*this->l_x_save)[i+k*(ngrid)]=xxsave[i];
191  for(int j=0;j<(*this->l_n_func);j++) {
192  (*this->l_y_save)[i+k*(ngrid)][j]=yysave[i][j];
193  }
194  }
195 
196  } else {
197 
198  /* iv_solver if we don't save */
199  this->oisp->solve_final_value
200  (xa,xb,h,(*this->l_n_func),y,y2,(*this->l_derivs));
201 
202  }
203 
204  /* Then we load values at the end of sub-interval */
205  if (k==(*this->l_mesh).size()-2) {
206  (*this->l_right_b)(xb,(*this->l_n_func),sx,y);
207  } else {
208  for(int i=0;i<(*this->l_n_func);i++) {
209  y[i]=sx[i+(*this->l_n_func)*(2+k)];
210  }
211  }
212 
213  /* Now we take the difference */
214  for(int i=0;i<(*this->l_n_func);i++) {
215  //sy[i+k*(*this->l_n_func)]=(y2[i]-y[i])/y[i];
216  sy[i+k*(*this->l_n_func)]= y2[i]-y[i];
217  }
218 
219  }
220 
221  /* Then load Extra boundary condition */
222  (*this->l_extra_b)(xb,(*this->l_n_func),sx,y);
223  for(int i=0;i<(*this->l_n_func);i++) {
224  sy[i+(int((*this->l_mesh).size()-1))*(*this->l_n_func)]=y[i];
225  }
226 
227  return 0;
228  }
229 
230 #endif
231 
232  };
233 
234 #ifndef DOXYGEN_NO_O2NS
235 }
236 #endif
237 
238 #endif
o2scl::ode_bv_multishoot::oisp
ode_iv_solve< func_t, vec_t, alloc_vec_t, alloc_t > * oisp
The initial value solver.
Definition: ode_bv_multishoot.h:121
boost::numeric::ublas::matrix< double >
o2scl::ode_iv_solve
Solve an initial-value ODE problems given an adaptive ODE stepper.
Definition: ode_iv_solve.h:103
o2scl::ode_iv_solve::solve_final_value
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, func_t &derivs)
Solve the initial-value problem to get the final value.
Definition: ode_iv_solve.h:197
o2scl::ode_bv_multishoot::l_y_save
mat_t * l_y_save
Desc.
Definition: ode_bv_multishoot.h:143
boost::numeric::ublas::vector< double >
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_bv_multishoot::save
bool save
Desc.
Definition: ode_bv_multishoot.h:145
o2scl::root
One-dimensional solver [abstract base].
Definition: root.h:48
o2scl::ode_bv_multishoot
Solve a ODE boundary value problem by multishooting.
Definition: ode_bv_multishoot.h:50
o2scl::ode_bv_multishoot::l_extra_b
func_t * l_extra_b
Desc.
Definition: ode_bv_multishoot.h:135
o2scl::ode_bv_multishoot::l_left_b
func_t * l_left_b
Desc.
Definition: ode_bv_multishoot.h:131
o2scl::ode_bv_multishoot::l_derivs
func_t * l_derivs
Desc.
Definition: ode_bv_multishoot.h:137
o2scl::mroot
Multidimensional root-finding [abstract base].
Definition: mroot.h:66
o2scl::ode_bv_multishoot::l_right_b
func_t * l_right_b
Desc.
Definition: ode_bv_multishoot.h:133
o2scl::ode_bv_multishoot::l_y_start
vec_t * l_y_start
Desc.
Definition: ode_bv_multishoot.h:129
o2scl::mm_funct
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef in src/base/mm_funct.h.
Definition: mm_funct.h:43
o2scl::ode_bv_multishoot::solve_fun
int solve_fun(size_t nv, const vec_t &sx, vec_t &sy)
Function to solve.
Definition: ode_bv_multishoot.h:148
o2scl::ode_bv_multishoot::l_x_save
vec_t * l_x_save
Desc.
Definition: ode_bv_multishoot.h:141
o2scl::ode_bv_multishoot::mrootp
gsl_mroot_hybrids< mm_funct<> > * mrootp
The equation solver.
Definition: ode_bv_multishoot.h:124
o2scl::ode_bv_multishoot::l_mesh
vec_t * l_mesh
Desc.
Definition: ode_bv_multishoot.h:127
o2scl::ode_bv_multishoot::l_n_func
int * l_n_func
Desc.
Definition: ode_bv_multishoot.h:139

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