TAPsODESolverRungeKutta4.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsODESolverRungeKutta4.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (07/10/2005)
UPDATE          (08/18/2010)
******************************************************************************/
#include "TAPsODESolverRungeKutta4.hpp"
// Using Inclusion Model (i.e. definitions are included in declarations)
//                       (this name.cpp is included in name.hpp)
// Each friend is defined directly inside its declaration.


BEGIN_NAMESPACE_TAPs__Simulation
//=============================================================================
//-----------------------------------------------------------------------------
template <typename T, typename T_SET>
ODESolverRungeKutta4<T, T_SET>::ODESolverRungeKutta4 () : ODESolver<T, T_SET>() {}
//-------------------------------------------------------------------
template <typename T, typename T_SET>
ODESolverRungeKutta4<T, T_SET>::~ODESolverRungeKutta4 ()    {}
//-------------------------------------------------------------------
template <typename T, typename T_SET>
void ODESolverRungeKutta4<T, T_SET>::SetSize ( int vectorSize )
{ 
    vK1.resize( vectorSize );
    vK2.resize( vectorSize );
    vK3.resize( vectorSize );
    vK4.resize( vectorSize );
}
//-------------------------------------------------------------------
template <typename T, typename T_SET>
void ODESolverRungeKutta4<T, T_SET>::Solve ( 
        VectorSet<T_SET> &x0,   // i/p: x(t_0) an initial state vector
        VectorSet<T_SET> &xEnd, // o/p: x(t_1) an end state vector
        T t0, T t1,             // i/p: starting and ending times
        DerivFunc dxdt,         // i/p: a derivative function
        void *userData          // o/p: array of user data
)
{
    bool restartNotRequired = true;
    T h = t1 - t0;
    T half = static_cast<T>( 0.5 );
    // Runge-Kutta of Order 4
    restartNotRequired = dxdt( t0, x0, vK1, userData );
    vK1 *= h;
    assert( restartNotRequired );
    restartNotRequired = dxdt( t0 + h*half, x0 + vK1*half, vK2, userData );
    vK2 *= h;
    assert( restartNotRequired );
    restartNotRequired = dxdt( t0 + h*half, x0 + vK2*half, vK3, userData );
    vK3 *= h;
    assert( restartNotRequired );
    restartNotRequired = dxdt( t0 + h, x0 + vK3, vK4, userData );
    vK4 *= h;
    assert( restartNotRequired );
    // Update the state
    xEnd = x0   + vK1/static_cast<T>(6) 
                + vK2/static_cast<T>(3)
                + vK3/static_cast<T>(3)
                + vK4/static_cast<T>(6);
}
//-------------------------------------------------------------------
//template <typename T, typename T_SET>
//void ODESolverRungeKutta4<T, T_SET>::Solve_ForPtrData ( 
//      VectorSetPtr<T_SET,T> &x0,      // i/p: x(t_0) an initial state vector
//      VectorSetPtr<T_SET,T> &xEnd,    // o/p: x(t_1) an end state vector
//      T t0, T t1,                     // i/p: starting and ending times
//      DerivFunc_ForPtrData dxdt,      // i/p: a derivative function
//      void *userData                  // o/p: array of user data
//)
//{
//}
//-----------------------------------------------------------------------------
//=============================================================================
END_NAMESPACE_TAPs__Simulation
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----