CppNoddy  0.90
Public Member Functions | List of all members
CppNoddy::ODE_IVP< _Type > Class Template Reference

A templated object for real/complex vector system of first-order ordinary differential equations. More...

#include <ODE_IVP.h>

Inheritance diagram for CppNoddy::ODE_IVP< _Type >:
CppNoddy::Uncopyable

Public Member Functions

 ODE_IVP (Equation< _Type > *equation_ptr, const double &x_init, const double &x_final, const std::size_t &num_of_steps)
 The class is defined by a vector function for the system. More...
 
 ~ODE_IVP ()
 
DenseVector< _Type > shoot4 (DenseVector< _Type > u)
 A fixed step 4th order Runge-Kutta method. More...
 
DenseVector< _Type > shoot45 (DenseVector< _Type > u, const double &tol, const double &h_init)
 A Runge-Kutta-Fehlberg integrator. More...
 
OneD_Node_Mesh< _Type > & get_mesh ()
 Return the history of the stepped solution. More...
 
unsigned & store_every ()
 Return a handle to the STORE_EVERY object. More...
 

Detailed Description

template<typename _Type>
class CppNoddy::ODE_IVP< _Type >

A templated object for real/complex vector system of first-order ordinary differential equations.

Definition at line 21 of file ODE_IVP.h.

Constructor & Destructor Documentation

template<typename _Type>
CppNoddy::ODE_IVP< _Type >::ODE_IVP ( Equation< _Type > *  equation_ptr,
const double &  x_init,
const double &  x_final,
const std::size_t &  num_of_steps 
)

The class is defined by a vector function for the system.

Parameters
equation_ptrA pointer to an inherited Equation object.
x_initThe starting point of the domain for the ODE.
x_finalThe end point of the domain for the ODE.
num_of_stepsA maximum/default number of steps to be taken.

Definition at line 20 of file ODE_IVP.cpp.

22  :
23  X_INIT( x1 ),
24  X_FINAL( x2 ),
25  H_INIT( ( x2 - x1 ) / num_of_points ),
26  N( num_of_points ),
27  p_EQUATION( ptr ),
28  STORE_EVERY( 1 )
29  {
30  p_EQUATION -> coord(0) = X_INIT;
31  }
template<typename _Type >
CppNoddy::ODE_IVP< _Type >::~ODE_IVP ( )

Definition at line 34 of file ODE_IVP.cpp.

35  {}

Member Function Documentation

template<typename _Type >
OneD_Node_Mesh< _Type > & CppNoddy::ODE_IVP< _Type >::get_mesh ( )

Return the history of the stepped solution.

Returns
A reference to the mesh solution

Definition at line 258 of file ODE_IVP.cpp.

Referenced by main().

259  {
260  return SOLN;
261  }
template<typename _Type>
DenseVector< _Type > CppNoddy::ODE_IVP< _Type >::shoot4 ( DenseVector< _Type >  u)

A fixed step 4th order Runge-Kutta method.

Parameters
uAn NVector of initial values.
Returns
An NVector of the final values.

Definition at line 38 of file ODE_IVP.cpp.

References h, CppNoddy::DenseVector< _Type >::push_back(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::set_nodes_vars(), CppNoddy::DenseVector< _Type >::size(), u, and CppNoddy::Example::z().

39  {
40 
41  double x = X_INIT;
42  const double h = H_INIT;
43  const double hby2 = h / 2.;
44  const double hby6 = h / 6.;
45  const double hby3 = h / 3.;
46  const int order = u.size();
47 
48  DenseVector<_Type> z( order, 0.0 ), k1( order, 0.0 ), k2( order, 0.0 ), k3( order, 0.0 ), k4( order, 0.0 );
49 
50  DenseVector< double > coords;
51  std::vector< DenseVector<_Type> > values;
52 
53  coords.push_back( x );
54  values.push_back( u );
55 
56  for ( unsigned i = 0; i < N; i++ )
57  {
58  // k1 = F(u,x)
59  p_EQUATION -> coord(0) = x;
60  p_EQUATION -> residual_fn( u, k1 );
61  z = u + k1 * hby2;
62 
63  x += hby2;
64 
65  // k2 = F(z,xhh)
66  p_EQUATION -> coord(0) = x;
67  p_EQUATION -> residual_fn( z, k2 );
68  z = u + k2 * hby2;
69 
70  // k3 = F(z,xhh)
71  p_EQUATION -> coord(0) = x;
72  p_EQUATION -> residual_fn( z, k3 );
73  z = u + k3 * h;
74 
75  x += hby2;
76 
77  // k4 = F(z,xh)
78  p_EQUATION -> coord(0) = x;
79  p_EQUATION -> residual_fn( z, k4 );
80  u += k1 * hby6 + k2 * hby3 + k3 * hby3 + k4 * hby6;
81 
82  if ( i % STORE_EVERY == 0 )
83  {
84  coords.push_back( x );
85  values.push_back( u );
86  }
87 
88  } //for loop stepping across domain
89 
90  // construct the solution mesh stored in this object
91  SOLN = OneD_Node_Mesh<_Type>( coords, p_EQUATION -> get_order() );
92  for ( unsigned i = 0; i < coords.size(); ++i )
93  {
94  // fill mesh
95  SOLN.set_nodes_vars( i, values[i] );
96  }
97  return u;
98  }
double z(const double &x)
Topography shape.
void set_nodes_vars(const std::size_t node, const DenseVector< _Type > &u)
Set the variables stored at A SPECIFIED node.
void push_back(const _Type &fill)
A pass-thru definition of push_back.
Definition: DenseVector.h:341
template<typename _Type>
DenseVector< _Type > CppNoddy::ODE_IVP< _Type >::shoot45 ( DenseVector< _Type >  u,
const double &  tol,
const double &  h_init 
)

A Runge-Kutta-Fehlberg integrator.

Parameters
uAn Nvector of initial values.
tolThe tolerance used in choosing the step length.
h_initThe initial step length.
Returns
An NVector of the final values.

Definition at line 101 of file ODE_IVP.cpp.

References h, CppNoddy::DenseVector< _Type >::push_back(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::set_nodes_vars(), CppNoddy::DenseVector< _Type >::size(), u, and CppNoddy::Example::z().

Referenced by main().

102  {
103  bool ok( false );
104  unsigned step = 0;
105  double x = X_INIT;
106  double h = h_init;
107  double c, diff;
108 
109  static const double X2 = 1. / 4.;
110  static const double X3 = 3. / 8.;
111  static const double X4 = 12. / 13.;
112  static const double X5 = 1.;
113  static const double X6 = 1. / 2.;
114 
115  static const double W21 = 1. / 4.;
116 
117  static const double W31 = 3. / 32.;
118  static const double W32 = 9. / 32.;
119 
120  static const double W41 = 1932. / 2197.;
121  static const double W42 = -7200. / 2197.;
122  static const double W43 = 7296. / 2197.;
123 
124  static const double W51 = 439. / 216.;
125  static const double W52 = -8.;
126  static const double W53 = 3680. / 513.;
127  static const double W54 = -845. / 4104;
128 
129  static const double W61 = -8. / 27.;
130  static const double W62 = 2.;
131  static const double W63 = -3544. / 2565.;
132  static const double W64 = 1859. / 4104.;
133  static const double W65 = -11. / 40.;
134 
135  static const double U1 = 25. / 216.;
136  static const double U3 = 1408. / 2565.;
137  static const double U4 = 2197. / 4104.;
138  static const double U5 = -1. / 5.;
139 
140  static const double Z1 = 16. / 135.;
141  static const double Z3 = 6656. / 12825.;
142  static const double Z4 = 28561. / 56430.;
143  static const double Z5 = -9. / 50.;
144  static const double Z6 = 2. / 55.;
145 
146  const unsigned order = u.size();
147 
148  DenseVector<_Type> z( order, 0.0 ), e( order, 0.0 ), k1( order, 0.0 ),
149  k2( order, 0.0 ), k3( order, 0.0 ), k4( order, 0.0 ), k5( order, 0.0 ), k6( order, 0.0 );
150 
151  DenseVector< double > coords;
152  std::vector< DenseVector<_Type> > values;
153 
154  coords.push_back( x );
155  values.push_back( u );
156 
157  do
158  {
159  step += 1;
160  // k1 = F(u,x)
161  p_EQUATION -> coord(0) = x;
162  p_EQUATION -> residual_fn( u, k1 );
163  k1 *= h;
164  z = u + k1 * W21;
165 
166  // k2 = F(z,x+X2*h)
167  p_EQUATION -> coord(0) = x + X2 * h;
168  p_EQUATION -> residual_fn( z, k2 );
169  k2 *= h;
170  z = u + k1 * W31 + k2 * W32;
171 
172  // k3 = F(z,x+X3*h)
173  p_EQUATION -> coord(0) = x + X3 * h;
174  p_EQUATION -> residual_fn( z, k3 );
175  k3 *= h;
176  z = u + k1 * W41 + k2 * W42 + k3 * W43;
177 
178  // k4 = F(z,x+X4*h)
179  p_EQUATION -> coord(0) = x + X4 * h;
180  p_EQUATION -> residual_fn( z, k4 );
181  k4 *= h;
182  z = u + k1 * W51 + k2 * W52 + k3 * W53 + k4 * W54;
183 
184  // k5 = F(z,x+X5*h)
185  p_EQUATION -> coord(0) = x + X5 * h;
186  p_EQUATION -> residual_fn( z, k5 );
187  k5 *= h;
188  z = u + k1 * W61 + k2 * W62 + k3 * W63 + k4 * W64 + k5 * W65;
189 
190  // k6 = F(z,x+X6*h)
191  p_EQUATION -> coord(0) = x + X6 * h;
192  p_EQUATION -> residual_fn( z, k6 );
193  k6 *= h;
194 
195  e = k1 * U1 + k3 * U3 + k4 * U4 + k5 * U5;
196  z = k1 * Z1 + k3 * Z3 + k4 * Z4 + k5 * Z5 + k6 * Z6;
197 
198  e -= z;
199 
200  // diff = ||e|| -- here use "abs" to deal with Complex systems
201  diff = e.inf_norm();
202 
203  c = sqrt( sqrt( tol * h / ( 2 * diff ) ) );
204  ok = true;
205 
206  // is the first step ok? or does it need reducing?
207 
208  if ( ( step == 1 ) && ( c < 1.0 ) )
209  {
210  // step needs reducing so start from initial value again
211  ok = false;
212  step = 1;
213  }
214 
215  if ( ok )
216  {
217  x += h;
218  u += z;
219 
220  if ( step % STORE_EVERY == 0 )
221  {
222  coords.push_back( x );
223  values.push_back( u );
224  }
225 
226  }
227 
228  h *= c;
229 
230  if ( x + h > X_FINAL )
231  {
232  h = ( X_FINAL - x );
233  }
234 
235  if ( step >= N )
236  {
237  std::string problem;
238  problem = "The ODE.shoot45 method reached the maximum \n";
239  problem += "number of steps specified by the user. \n";
240  throw ExceptionRuntime( problem );
241  }
242 
243  }
244  while ( std::abs( x - X_FINAL ) > tol ); // end loop stepping across domain
245 
246  // construct the solution mesh stored in this object
247  SOLN = OneD_Node_Mesh<_Type>( coords, p_EQUATION -> get_order() );
248  for ( unsigned i = 0; i < coords.size(); ++i )
249  {
250  // fill mesh
251  SOLN.set_nodes_vars( i, values[i] );
252  }
253 
254  return u;
255  }
double z(const double &x)
Topography shape.
void set_nodes_vars(const std::size_t node, const DenseVector< _Type > &u)
Set the variables stored at A SPECIFIED node.
void push_back(const _Type &fill)
A pass-thru definition of push_back.
Definition: DenseVector.h:341
template<typename _Type >
unsigned & CppNoddy::ODE_IVP< _Type >::store_every ( )

Return a handle to the STORE_EVERY object.

Returns
A reference to the STORE_EVERY variable

Definition at line 264 of file ODE_IVP.cpp.

Referenced by main().

265  {
266  return STORE_EVERY;
267  }

The documentation for this class was generated from the following files:

© 2012

R.E. Hewitt