CppNoddy  0.85
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Groups Pages
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.
 ~ODE_IVP ()
DenseVector< _Type > shoot4 (DenseVector< _Type > u)
 A fixed step 4th order Runge-Kutta method.
DenseVector< _Type > shoot45 (DenseVector< _Type > u, const double &tol, const double &h_init)
 A Runge-Kutta-Fehlberg integrator.
OneD_Node_Mesh< _Type > & get_mesh ()
 Return the history of the stepped solution.
unsigned & store_every ()
 Return a handle to the STORE_EVERY object.

Additional Inherited Members

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.

:
X_INIT( x1 ),
X_FINAL( x2 ),
H_INIT( ( x2 - x1 ) / num_of_points ),
N( num_of_points ),
p_EQUATION( ptr ),
STORE_EVERY( 1 )
{
p_EQUATION -> coord(0) = X_INIT;
}
template<typename _Type >
CppNoddy::ODE_IVP< _Type >::~ODE_IVP ( )

Definition at line 34 of file ODE_IVP.cpp.

{}

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 257 of file ODE_IVP.cpp.

{
return SOLN;
}
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().

{
double x = X_INIT;
const double h = H_INIT;
const double hby2 = h / 2.;
const double hby6 = h / 6.;
const double hby3 = h / 3.;
const int order = u.size();
DenseVector<_Type> z( order, 0.0 ), k1( order, 0.0 ), k2( order, 0.0 ), k3( order, 0.0 ), k4( order, 0.0 );
DenseVector< double > coords;
std::vector< DenseVector<_Type> > values;
coords.push_back( x );
values.push_back( u );
for ( unsigned i = 0; i < N; i++ )
{
// k1 = F(u,x)
p_EQUATION -> coord(0) = x;
p_EQUATION -> residual_fn( u, k1 );
z = u + k1 * hby2;
x += hby2;
// k2 = F(z,xhh)
p_EQUATION -> coord(0) = x;
p_EQUATION -> residual_fn( z, k2 );
z = u + k2 * hby2;
// k3 = F(z,xhh)
p_EQUATION -> coord(0) = x;
p_EQUATION -> residual_fn( z, k3 );
z = u + k3 * h;
x += hby2;
// k4 = F(z,xh)
p_EQUATION -> coord(0) = x;
p_EQUATION -> residual_fn( z, k4 );
u += k1 * hby6 + k2 * hby3 + k3 * hby3 + k4 * hby6;
if ( i % STORE_EVERY == 0 )
{
coords.push_back( x );
values.push_back( u );
}
} //for loop stepping across domain
// construct the solution mesh stored in this object
SOLN = OneD_Node_Mesh<_Type>( coords, p_EQUATION -> get_order() );
for ( unsigned i = 0; i < coords.size(); ++i )
{
// fill mesh
SOLN.set_nodes_vars( i, values[i] );
}
return u;
}
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 100 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().

{
bool ok( false );
unsigned step = 0;
double x = X_INIT;
double h = h_init;
double c, diff;
static const double X2 = 1. / 4.;
static const double X3 = 3. / 8.;
static const double X4 = 12. / 13.;
static const double X5 = 1.;
static const double X6 = 1. / 2.;
static const double W21 = 1. / 4.;
static const double W31 = 3. / 32.;
static const double W32 = 9. / 32.;
static const double W41 = 1932. / 2197.;
static const double W42 = -7200. / 2197.;
static const double W43 = 7296. / 2197.;
static const double W51 = 439. / 216.;
static const double W52 = -8.;
static const double W53 = 3680. / 513.;
static const double W54 = -845. / 4104;
static const double W61 = -8. / 27.;
static const double W62 = 2.;
static const double W63 = -3544. / 2565.;
static const double W64 = 1859. / 4104.;
static const double W65 = -11. / 40.;
static const double U1 = 25. / 216.;
static const double U3 = 1408. / 2565.;
static const double U4 = 2197. / 4104.;
static const double U5 = -1. / 5.;
static const double Z1 = 16. / 135.;
static const double Z3 = 6656. / 12825.;
static const double Z4 = 28561. / 56430.;
static const double Z5 = -9. / 50.;
static const double Z6 = 2. / 55.;
const unsigned order = u.size();
DenseVector<_Type> z( order, 0.0 ), e( order, 0.0 ), k1( order, 0.0 ),
k2( order, 0.0 ), k3( order, 0.0 ), k4( order, 0.0 ), k5( order, 0.0 ), k6( order, 0.0 );
DenseVector< double > coords;
std::vector< DenseVector<_Type> > values;
coords.push_back( x );
values.push_back( u );
do
{
step += 1;
// k1 = F(u,x)
p_EQUATION -> coord(0) = x;
p_EQUATION -> residual_fn( u, k1 );
k1 *= h;
z = u + k1 * W21;
// k2 = F(z,x+X2*h)
p_EQUATION -> coord(0) = x + X2 * h;
p_EQUATION -> residual_fn( z, k2 );
k2 *= h;
z = u + k1 * W31 + k2 * W32;
// k3 = F(z,x+X3*h)
p_EQUATION -> coord(0) = x + X3 * h;
p_EQUATION -> residual_fn( z, k3 );
k3 *= h;
z = u + k1 * W41 + k2 * W42 + k3 * W43;
// k4 = F(z,x+X4*h)
p_EQUATION -> coord(0) = x + X4 * h;
p_EQUATION -> residual_fn( z, k4 );
k4 *= h;
z = u + k1 * W51 + k2 * W52 + k3 * W53 + k4 * W54;
// k5 = F(z,x+X5*h)
p_EQUATION -> coord(0) = x + X5 * h;
p_EQUATION -> residual_fn( z, k5 );
k5 *= h;
z = u + k1 * W61 + k2 * W62 + k3 * W63 + k4 * W64 + k5 * W65;
// k6 = F(z,x+X6*h)
p_EQUATION -> coord(0) = x + X6 * h;
p_EQUATION -> residual_fn( z, k6 );
k6 *= h;
e = k1 * U1 + k3 * U3 + k4 * U4 + k5 * U5;
z = k1 * Z1 + k3 * Z3 + k4 * Z4 + k5 * Z5 + k6 * Z6;
e -= z;
// diff = ||e|| -- here use "abs" to deal with Complex systems
diff = e.inf_norm();
c = sqrt( sqrt( tol * h / ( 2 * diff ) ) );
ok = true;
// is the first step ok? or does it need reducing?
if ( ( step == 1 ) && ( c < 1.0 ) )
{
// step needs reducing so start from initial value again
ok = false;
step = 1;
}
if ( ok )
{
x += h;
u += z;
if ( step % STORE_EVERY == 0 )
{
coords.push_back( x );
values.push_back( u );
}
}
h *= c;
if ( x + h > X_FINAL )
{
h = ( X_FINAL - x );
}
if ( step >= N )
{
std::string problem;
problem = "The ODE.shoot45 method reached the maximum \n";
problem += "number of steps specified by the user. \n";
throw ExceptionRuntime( problem );
}
}
while ( std::abs( x - X_FINAL ) > tol ); // end loop stepping across domain
// construct the solution mesh stored in this object
SOLN = OneD_Node_Mesh<_Type>( coords, p_EQUATION -> get_order() );
for ( unsigned i = 0; i < coords.size(); ++i )
{
// fill mesh
SOLN.set_nodes_vars( i, values[i] );
}
return u;
}
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 263 of file ODE_IVP.cpp.

{
return STORE_EVERY;
}

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

© 2012

R.E. Hewitt