ode_iv_table.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_IV_TABLE_H
24 #define O2SCL_ODE_IV_TABLE_H
25 
26 /** \file ode_iv_table.h
27  \brief File defining \ref o2scl::ode_iv_table
28 */
29 
30 #include <o2scl/astep.h>
31 #include <o2scl/astep_gsl.h>
32 #include <o2scl/ode_iv_solve.h>
33 
34 #ifndef DOXYGEN_NO_O2NS
35 namespace o2scl {
36 #endif
37 
38  /** \brief Solve an initial-value ODE problem and store
39  the result in a \ref table object
40 
41  This class is experimental.
42 
43  \future It would be nice not to have to copy the results from a
44  matrix into a table, but this may require a nontrivial
45  modification of the ODE solvers and/or the table class.
46 
47  \future One possible idea is to redo the name specification as
48  a separate function, which allows one to either specify
49  prefixes or full column names. We also need to figure out
50  how to handle column units.
51  */
52  template<class func_t=ode_funct<>,
53  class vec_t=ubvector, class alloc_vec_t=ubvector,
54  class alloc_t=ubvector_alloc> class ode_iv_table :
55  public ode_iv_solve<func_t,vec_t,alloc_vec_t,alloc_t> {
56 
57  public:
58 
59  /// Desc
60  int solve_grid_table(size_t n, vec_t &ystart, table<> &t, std::string x_col,
61  std::string y_prefix, std::string dydx_prefix,
62  std::string yerr_prefix, func_t &derivs) {
63  std::string cname;
64 
65  size_t n_sol=t.get_nlines();
66  double x0=t.get(x_col,0);
67  double x1=t.get(x_col,n_sol-1);
68  double h=t.get(x_col,1)-x0;
69 
70  // Allocate vectors and matrices for solve_grid()
71  alloc_vec_t x_sol;
72  this->ao.allocate(x_sol,n_sol);
73  ubmatrix ysol(n_sol,n), dydx_sol(n_sol,n), yerr_sol(n_sol,n);
74 
75  // Copy over x_grid from table
76  ubvector &ob=t.get_column(x_col);
77  vector_copy(n_sol,x_col,x_sol);
78 
79  // Perform solution
80  this->template solve_grid<ubmatrix,ubmatrix_row>
81  (x0,x1,h,n,ystart,n_sol,x_sol,ysol,dydx_sol,yerr_sol,derivs);
82 
83  // Pointers to columns in table
84  std::vector<ubvector *> yt, dyt, errt;
85 
86  // Create new columns if necessary and get pointers
87  for(size_t i=0;i<n;i++) {
88 
89  // TABLE FIXME
90 
91  /*
92  cname=y_prefix+itos(i);
93  if (!t.is_column(cname)) t.new_column(cname);
94  yt.push_back(&t.get_column(cname));
95 
96  cname=dydx_prefix+itos(i);
97  if (!t.is_column(cname)) t.new_column(cname);
98  dyt.push_back(&t.get_column(cname));
99 
100  cname=yerr_prefix+itos(i);
101  if (!t.is_column(cname)) t.new_column(cname);
102  errt.push_back(&t.get_column(cname));
103 
104  // Now copy data over
105  for(size_t j=0;j<n_sol;j++) {
106  (*yt[i])[j]=ysol[j][i];
107  (*dyt[i])[j]=dydx_sol[j][i];
108  (*errt[i])[j]=yerr_sol[j][i];
109  }
110  */
111 
112  }
113 
114  return 0;
115  }
116 
117  /// Desc
118  int solve_store_table(double x0, double x1, double h, size_t n,
119  vec_t &ystart, size_t &n_sol, table<> &t,
120  std::string x_col, std::string y_prefix,
121  std::string dydx_prefix, std::string yerr_prefix,
122  func_t &derivs) {
123 
124  std::string cname;
125  if (t.get_nlines()<n_sol) t.set_nlines(n_sol);
126 
127  // Allocate vectors and matrices for solve_store()
128  alloc_vec_t x_sol;
129  this->ao.allocate(x_sol,n_sol);
130  ubmatrix ysol(n_sol,n), dydx_sol(n_sol,n), yerr_sol(n_sol,n);
131 
132  // Perform solution
133  this->template solve_store<ubmatrix,ubmatrix_row>
134  (x0,x1,h,n,ystart,n_sol,x_sol,ysol,dydx_sol,yerr_sol,derivs);
135 
136  // Pointers to columns in table
137  std::vector<ubvector *> yt, dyt, errt;
138 
139  if (!t.is_column(x_col)) t.new_column(x_col);
140  ubvector &x_vec=t.get_column(x_col);
141 
142  // Create new columns if necessary and get pointers
143  for(size_t i=0;i<n;i++) {
144 
145  // TABLE FIXME
146 
147  /*
148 
149  cname=y_prefix+itos(i);
150  if (!t.is_column(cname)) t.new_column(cname);
151  yt.push_back(&t.get_column(cname));
152 
153  cname=dydx_prefix+itos(i);
154  if (!t.is_column(cname)) t.new_column(cname);
155  dyt.push_back(&t.get_column(cname));
156 
157  cname=yerr_prefix+itos(i);
158  if (!t.is_column(cname)) t.new_column(cname);
159  errt.push_back(&t.get_column(cname));
160 
161  // Now copy data over
162  for(size_t j=0;j<n_sol;j++) {
163  if (i==0) x_vec[j]=x_sol[j];
164  (*yt[i])[j]=ysol[j][i];
165  (*dyt[i])[j]=dydx_sol[j][i];
166  (*errt[i])[j]=yerr_sol[j][i];
167  }
168 
169  */
170 
171  }
172 
173  return 0;
174  }
175 
176  };
177 
178 #ifndef DOXYGEN_NO_O2NS
179 }
180 #endif
181 
182 #endif
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
Definition: table.h:414
Solve an initial-value ODE problems given an adaptive ODE stepper.
Definition: ode_iv_solve.h:91
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 get_nlines() const
Return the number of lines.
Definition: table.h:457
Data table table class.
Definition: table.h:49
void set_nlines(size_t il)
Set the number of lines.
Definition: table.h:467
const vec_t & get_column(std::string scol) const
Returns a reference to the column named col. .
Definition: table.h:619
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
int solve_store_table(double x0, double x1, double h, size_t n, vec_t &ystart, size_t &n_sol, table<> &t, std::string x_col, std::string y_prefix, std::string dydx_prefix, std::string yerr_prefix, func_t &derivs)
Desc.
Definition: ode_iv_table.h:118
Solve an initial-value ODE problem and store the result in a table object.
Definition: ode_iv_table.h:54
int solve_grid_table(size_t n, vec_t &ystart, table<> &t, std::string x_col, std::string y_prefix, std::string dydx_prefix, std::string yerr_prefix, func_t &derivs)
Desc.
Definition: ode_iv_table.h:60
void new_column(std::string head)
Add a new column owned by the table table .
Definition: table.h:690
bool is_column(std::string scol) const
Return true if scol is a column in the current table table .
Definition: table.h:899
static const double x1[5]
Definition: inte_qng_gsl.h:48

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