ode_iv_table.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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,
61  std::string x_col, std::string y_prefix,
62  std::string dydx_prefix, std::string yerr_prefix,
63  func_t &derivs) {
64 
65  std::string cname;
66 
67  size_t n_sol=t.get_nlines();
68  double x0=t.get(x_col,0);
69  double x1=t.get(x_col,n_sol-1);
70  double h=t.get(x_col,1)-x0;
71 
72  // Allocate vectors and matrices for solve_grid()
73  alloc_vec_t x_sol;
74  this->ao.allocate(x_sol,n_sol);
75  ubmatrix ysol(n_sol,n), dydx_sol(n_sol,n), yerr_sol(n_sol,n);
76 
77  // Copy over x_grid from table
78  ubvector &ob=t.get_column(x_col);
79  vector_copy(n_sol,x_col,x_sol);
80 
81  // Perform solution
82  this->template solve_grid<ubmatrix,ubmatrix_row>
83  (x0,x1,h,n,ystart,n_sol,x_sol,ysol,dydx_sol,yerr_sol,derivs);
84 
85  // Pointers to columns in table
86  std::vector<ubvector *> yt, dyt, errt;
87 
88  // Create new columns if necessary and get pointers
89  for(size_t i=0;i<n;i++) {
90 
91  // TABLE FIXME
92 
93  /*
94  cname=y_prefix+itos(i);
95  if (!t.is_column(cname)) t.new_column(cname);
96  yt.push_back(&t.get_column(cname));
97 
98  cname=dydx_prefix+itos(i);
99  if (!t.is_column(cname)) t.new_column(cname);
100  dyt.push_back(&t.get_column(cname));
101 
102  cname=yerr_prefix+itos(i);
103  if (!t.is_column(cname)) t.new_column(cname);
104  errt.push_back(&t.get_column(cname));
105 
106  // Now copy data over
107  for(size_t j=0;j<n_sol;j++) {
108  (*yt[i])[j]=ysol[j][i];
109  (*dyt[i])[j]=dydx_sol[j][i];
110  (*errt[i])[j]=yerr_sol[j][i];
111  }
112  */
113 
114  }
115 
116  return 0;
117  }
118 
119  /// Desc
120  int solve_store_table(double x0, double x1, double h, size_t n,
121  vec_t &ystart, size_t &n_sol, table<> &t,
122  std::string x_col, std::string y_prefix,
123  std::string dydx_prefix, std::string yerr_prefix,
124  func_t &derivs) {
125 
126  std::string cname;
127  if (t.get_nlines()<n_sol) t.set_nlines(n_sol);
128 
129  // Allocate vectors and matrices for solve_store()
130  alloc_vec_t x_sol;
131  this->ao.allocate(x_sol,n_sol);
132  ubmatrix ysol(n_sol,n), dydx_sol(n_sol,n), yerr_sol(n_sol,n);
133 
134  // Perform solution
135  this->template solve_store<ubmatrix,ubmatrix_row>
136  (x0,x1,h,n,ystart,n_sol,x_sol,ysol,dydx_sol,yerr_sol,derivs);
137 
138  // Pointers to columns in table
139  std::vector<ubvector *> yt, dyt, errt;
140 
141  if (!t.is_column(x_col)) t.new_column(x_col);
142  ubvector &x_vec=t.get_column(x_col);
143 
144  // Create new columns if necessary and get pointers
145  for(size_t i=0;i<n;i++) {
146 
147  // TABLE FIXME
148 
149  /*
150 
151  cname=y_prefix+itos(i);
152  if (!t.is_column(cname)) t.new_column(cname);
153  yt.push_back(&t.get_column(cname));
154 
155  cname=dydx_prefix+itos(i);
156  if (!t.is_column(cname)) t.new_column(cname);
157  dyt.push_back(&t.get_column(cname));
158 
159  cname=yerr_prefix+itos(i);
160  if (!t.is_column(cname)) t.new_column(cname);
161  errt.push_back(&t.get_column(cname));
162 
163  // Now copy data over
164  for(size_t j=0;j<n_sol;j++) {
165  if (i==0) x_vec[j]=x_sol[j];
166  (*yt[i])[j]=ysol[j][i];
167  (*dyt[i])[j]=dydx_sol[j][i];
168  (*errt[i])[j]=yerr_sol[j][i];
169  }
170 
171  */
172 
173  }
174 
175  return 0;
176  }
177 
178  };
179 
180 #ifndef DOXYGEN_NO_O2NS
181 }
182 #endif
183 
184 #endif
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
Definition: table.h:422
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:465
Data table table class.
Definition: table.h:49
void set_nlines(size_t il)
Set the number of lines.
Definition: table.h:475
const vec_t & get_column(std::string scol) const
Returns a reference to the column named col. .
Definition: table.h:671
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:120
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:742
bool is_column(std::string scol) const
Return true if scol is a column in the current table table .
Definition: table.h:953
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).