ode_bv_multishoot.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2008, 2009, 2010, 2011, Julien Garaud and
5  Andrew W. Steiner
6 
7  This file is part of O2scl.
8 
9  O2scl is free software; you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version.
13 
14  O2scl is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
21 
22  -------------------------------------------------------------------
23 */
24 #ifndef O2SCL_ODE_BV_MULTISHOOT_H
25 #define O2SCL_ODE_BV_MULTISHOOT_H
26 
27 /** \file ode_bv_multishoot.h
28  \brief File defining \ref o2scl::ode_bv_multishoot
29 */
30 
31 #include <string>
32 #include <o2scl/astep.h>
33 #include <o2scl/astep_gsl.h>
34 #include <o2scl/ode_iv_solve.h>
35 #include <o2scl/gsl_mroot_hybrids.h>
36 
37 #ifndef DOXYGEN_NO_O2NS
38 namespace o2scl {
39 #endif
40 
41  /** \brief Solve a ODE boundary value problem by multishooting
42 
43  This class is experimental.
44 
45  \todo Improve documentation a little and create
46  testing code
47  */
48  template <class func_t=ode_funct<>,
49  class vec_t=ubvector, class alloc_vec_t=ubvector,
50  class alloc_t=ubvector_alloc,class vec_int_t=ubvector_int_base,
51  class mat_t=ubmatrix> class ode_bv_multishoot {
52 
53  public:
54 
56  oisp=&def_ois;
57  mrootp=&def_mroot;
58  }
59 
60  virtual ~ode_bv_multishoot() {
61  }
62 
63  /* Solve the boundary value problem on a mesh */
64  virtual int solve(vec_t &mesh, int &n_func, vec_t &y_start,
65  func_t &left_b, func_t &right_b, func_t &extra_b,
66  func_t &derivs, vec_t &x_save, mat_t &y_save) {
67 
68  /* Make copies of the input for later access */
69  this->l_mesh=&mesh;
70  this->l_n_func=&n_func;
71  this->l_y_start=&y_start;
72  this->l_left_b=&left_b;
73  this->l_right_b=&right_b;
74  this->l_extra_b=&extra_b;
75  this->l_derivs=&derivs;
76  this->l_x_save=&x_save;
77  this->l_y_save=&y_save;
78 
79  /* vector of variables */
80  int nsolve=y_start.size();
81  ubvector sx(nsolve),sy(nsolve);
82  sx=y_start;
83 
84  /* Equation solver */
85  mm_funct_mfptr<ode_bv_multishoot<func_t,vec_t,
86  alloc_vec_t,alloc_t,vec_int_t> >
87  mfm(this,&ode_bv_multishoot<func_t,vec_t,
88  alloc_vec_t,alloc_t,vec_int_t>::solve_fun);
89 
90  /* Run multishooting and save at the last step */
91  this->save=false;
92  int ret=this->mrootp->msolve(nsolve,sx,mfm);
93  this->save=true;
94  solve_fun(nsolve,sx,sy);
95 
96  return ret;
97  }
98 
99  /* Set initial value solver */
101  oisp=&ois;
102  return 0;
103  }
104 
105  /* Set equation solver */
106  int set_mroot(mroot<mm_funct<> > &root) {
107  mrootp=&root;
108  return 0;
109  }
110 
111  /* Default initial value solver */
113 
114  /* Default equation solver */
115  gsl_mroot_hybrids<mm_funct<> > def_mroot;
116 
117 #ifndef DOXYGEN_INTERNAL
118 
119  protected:
120 
121  /// The initial value solver
123 
124  /// The equation solver
125  gsl_mroot_hybrids<mm_funct<> > *mrootp;
126 
127  /// Desc
128  vec_t *l_mesh;
129  /// Desc
130  vec_t *l_y_start;
131  /// Desc
132  func_t *l_left_b;
133  /// Desc
134  func_t *l_right_b;
135  /// Desc
136  func_t *l_extra_b;
137  /// Desc
138  func_t *l_derivs;
139  /// Desc
140  int *l_n_func;
141  /// Desc
142  vec_t *l_x_save;
143  /// Desc
144  mat_t *l_y_save;
145  /// Desc
146  bool save;
147 
148  /// Function to solve
149  int solve_fun(size_t nv, const vec_t &sx, vec_t &sy) {
150 
151  double xa,xb=0.0,h;
152  ubvector y((*this->l_n_func)),y2((*this->l_n_func));
153 
154  /* We update y_start in order that derivs knows
155  all the values of parameters */
156  for(size_t i=0;i<(*this->l_y_start).size();i++) {
157  (*this->l_y_start)[i]=sx[i];
158  }
159 
160  /* A loop on each subinterval */
161  for(size_t k=0;k<(*this->l_mesh).size()-1;k++) {
162 
163  xa=(*this->l_mesh)[k];
164  xb=(*this->l_mesh)[k+1];
165  h=(xb-xa)/100.0;
166 
167  /* We load function's value at the left point of the sub-interval */
168  if (k==0) {
169  (*this->l_left_b)(xa,(*this->l_n_func),sx,y);
170  } else {
171  for(int i=0;i<(*this->l_n_func);i++)
172  y[i]=sx[i+(*this->l_n_func)*(1+k)];
173  }
174 
175  if (this->save) {
176 
177  /* iv_solver if we save */
178  int ngrid=((*this->l_x_save).size()-1)/((*this->l_mesh).size()-1)+1;
179  ubvector xxsave(ngrid);
180  ubmatrix yysave(ngrid,(*this->l_n_func));
181 
182  if (k!=((*this->l_mesh).size()-2)) {
183  xb=(*this->l_mesh)[k+1]-((*this->l_mesh)[k+1]-
184  (*this->l_mesh)[k])/ngrid;
185  }
186 
187  this->oisp->solve_grid(xa,xb,h,(*this->l_n_func),y,ngrid,
188  xxsave,yysave,(*this->l_derivs));
189 
190  for(int i=0;i<ngrid;i++) {
191  (*this->l_x_save)[i+k*(ngrid)]=xxsave[i];
192  for(int j=0;j<(*this->l_n_func);j++) {
193  (*this->l_y_save)[i+k*(ngrid)][j]=yysave[i][j];
194  }
195  }
196 
197  } else {
198 
199  /* iv_solver if we don't save */
200  this->oisp->solve_final_value
201  (xa,xb,h,(*this->l_n_func),y,y2,(*this->l_derivs));
202 
203  }
204 
205  /* Then we load values at the end of sub-interval */
206  if (k==(*this->l_mesh).size()-2) {
207  (*this->l_right_b)(xb,(*this->l_n_func),sx,y);
208  } else {
209  for(int i=0;i<(*this->l_n_func);i++) {
210  y[i]=sx[i+(*this->l_n_func)*(2+k)];
211  }
212  }
213 
214  /* Now we take the difference */
215  for(int i=0;i<(*this->l_n_func);i++) {
216  //sy[i+k*(*this->l_n_func)]=(y2[i]-y[i])/y[i];
217  sy[i+k*(*this->l_n_func)]= y2[i]-y[i];
218  }
219 
220  }
221 
222  /* Then load Extra boundary condition */
223  (*this->l_extra_b)(xb,(*this->l_n_func),sx,y);
224  for(int i=0;i<(*this->l_n_func);i++) {
225  sy[i+(int((*this->l_mesh).size()-1))*(*this->l_n_func)]=y[i];
226  }
227 
228  return 0;
229  }
230 
231 #endif
232 
233  };
234 
235 #ifndef DOXYGEN_NO_O2NS
236 }
237 #endif
238 
239 #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_iv_solve< func_t, vec_t, alloc_vec_t, alloc_t > * oisp
The initial value solver.
Solve a ODE boundary value problem by multishooting.
gsl_mroot_hybrids< mm_funct<> > * mrootp
The equation solver.
Multidimensional root-finding [abstract base].
Definition: mroot.h:66
One-dimensional solver [abstract base].
Definition: root.h:47
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:181
int solve_fun(size_t nv, const vec_t &sx, vec_t &sy)
Function to solve.
int solve_grid(double h, size_t n, size_t nsol, vec_t &xsol, mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol, func_t &derivs)
Solve the initial-value problem from x0 to x1 over a grid storing derivatives and errors...
Definition: ode_iv_solve.h:612

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