23 #ifndef O2SCL_ODE_IT_SOLVE_H 24 #define O2SCL_ODE_IT_SOLVE_H 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> 37 #include <o2scl/misc.h> 38 #include <o2scl/test_mgr.h> 39 #include <o2scl/linear_solver.h> 40 #include <o2scl/vector.h> 42 #ifndef DOXYGEN_NO_O2NS 47 typedef std::function<double
67 class solver_mat_t=boost::numeric::ublas::matrix<double> >
75 eps_rel=sqrt(std::numeric_limits<double>::epsilon());
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) {
145 typedef std::function<double
146 (
size_t,
size_t,
double,matrix_row_t &)> ode_it_dfunct;
148 ode_it_dfunct d2_derivs=std::bind
149 (std::mem_fn<
double(
size_t,
size_t,
double,matrix_row_t &)>
151 std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
152 ode_it_dfunct d2_left=std::bind
153 (std::mem_fn<
double(
size_t,
size_t,
double,matrix_row_t &)>
155 std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
156 ode_it_dfunct d2_right=std::bind
157 (std::mem_fn<
double(
size_t,
size_t,
double,matrix_row_t &)>
159 std::placeholders::_2,std::placeholders::_3,std::placeholders::_4);
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);
183 template<
class dfunc_t>
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) {
193 size_t nb_right=n_eq-nb_left;
196 size_t nvars=n_grid*n_eq;
199 for(
size_t it=0;done==
false && it<
niter;it++) {
203 for(
size_t i=0;i<nvars;i++) {
204 for(
size_t j=0;j<nvars;j++) {
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);
222 for(
size_t k=0;k<n_grid-1;k++) {
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);
229 for(
size_t i=0;i<n_eq;i++) {
231 rhs[ix]=y(k,i)-y(kp1,i)+(x[kp1]-x[k])*
232 (derivs(i,tx,ykp1)+derivs(i,tx,yk))/2.0;
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;
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;
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);
256 rhs[ix]=-right(i,x[n_grid-1],ylast);
258 for(
size_t j=0;j<n_eq;j++) {
259 mat(ix,lhs+j)=d_right(i,j,x[n_grid-1],ylast);
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) <<
" ";
274 std::cout << std::endl;
276 std::cout <<
"Deviations:" << std::endl;
277 for(
size_t i=0;i<nvars;i++) {
278 std::cout << rhs[i] << std::endl;
282 if (make_mats)
return 0;
287 std::cout <<
"Corrections:" << std::endl;
288 for(
size_t i=0;i<nvars;i++) {
289 std::cout << dy[i] << std::endl;
291 std::cout <<
"Press a key and press enter to continue: " 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];
313 std::cout <<
"ode_it_solve: " << it <<
" " << std::sqrt(res) <<
" " 314 << tol_rel << std::endl;
317 std::cout <<
"Press a key and type enter to continue. ";
323 if (std::sqrt(res)<=tol_rel) done=
true;
327 O2SCL_ERR(
"Exceeded number of iterations in solve().",
341 func_t *fl, *fr, *fd;
352 virtual double fd_left(
size_t ieq,
size_t ivar,
double x, matrix_row_t &y) {
356 double h=eps_rel*fabs(y[ivar]);
375 virtual double fd_right(
size_t ieq,
size_t ivar,
double x, matrix_row_t &y) {
379 double h=eps_rel*fabs(y[ivar]);
399 virtual double fd_derivs(
size_t ieq,
size_t ivar,
double x,
404 double h=eps_rel*fabs(y[ivar]);
421 #ifndef DOXYGEN_NO_O2NS
std::function< double(size_t, double, boost::numeric::ublas::matrix_row< boost::numeric::ublas::matrix< double > > &)> ode_it_funct
Function for iterative solving of ODEs.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
ODE solver using a generic linear solver to solve finite-difference equations.
int set_solver(o2scl_linalg::linear_solver< solver_vec_t, solver_mat_t > &ls)
Set the linear solver.
double alpha
Size of correction to apply (default 1.0)
exceeded max number of iterations
int verbose
Set level of output (default 0)
o2scl_linalg::linear_solver< solver_vec_t, solver_mat_t > * solver
Solver.
virtual double fd_left(size_t ieq, size_t ivar, double x, matrix_row_t &y)
Compute the derivatives of the LHS boundary conditions.
double eps_rel
Stepsize for finite differencing (default )
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.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
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.
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
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.
virtual void solve(size_t n, mat_t &a, vec_t &b, vec_t &x)=0
Solve square linear system of size n.
double tol_rel
Tolerance (default )
double eps_min
Minimum stepsize for finite differencing (default )
virtual double fd_right(size_t ieq, size_t ivar, double x, matrix_row_t &y)
Compute the derivatives of the RHS boundary conditions.
o2scl_linalg::linear_solver_HH< solver_vec_t, solver_mat_t > def_solver
Default linear solver.
size_t niter
Maximum number of iterations (default 30)