Back to index

kdeartwork  4.3.2
Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes
RkOdeSolver< T, D > Class Template Reference

Solver class to integrate a first-order ordinary differential equation (ODE) by means of a 6. More...

#include <rkodesolver.h>

Inheritance diagram for RkOdeSolver< T, D >:
Inheritance graph
[legend]
Collaboration diagram for RkOdeSolver< T, D >:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 RkOdeSolver (const T &x, const Matrix< T, D, 1 > &y, const T &dx=0.0, const T &eps=1e-6)
 Constructor.
virtual ~RkOdeSolver (void)
 Destructor.
void integrate (const T &dx)
 Integrates the ordinary differential equation from the current x value to x+dx.
const T & X (void) const
 Get current x value.
void X (const T &a)
 Set current x value.
const Matrix< T, D, 1 > & Y (void) const
 Get current y value.
void Y (const Matrix< T, D, 1 > &a)
 Set current y values.
const Matrix< T, D, 1 > & dYdX (void) const
 Get current dy/dx value.
const T & dX (void) const
 Get current estimated step size dX.
void dX (const T &a)
 Set estimated step size dX.
const T & Eps (void) const
 Get current presision.
void Eps (const T &a)
 Set estimated presision.

Protected Member Functions

virtual Matrix< T, D, 1 > f (const T &x, const Matrix< T, D, 1 > &y) const =0
 ODE function.

Private Member Functions

bool rkStepCheck (const T &dx)
 Perform one integration step with a tolerable relative error given by ::mErr.
Matrix< T, D, 1 > rkStep (const T &dx, Matrix< T, D, 1 > &yerr) const
 Perform one Runge-Kutta integration step forward with step size ::m_step.

Private Attributes

m_x
 current x value
Matrix< T, D, 1 > m_y
 current y value
Matrix< T, D, 1 > m_dydx
 current value of dy/dx
m_eps
 allowed relative error
m_step
 estimated step size for next Runge-Kutta step

Detailed Description

template<typename T, int D>
class RkOdeSolver< T, D >

Solver class to integrate a first-order ordinary differential equation (ODE) by means of a 6.

order Runge-Kutta method.

See the article about the Cash-Karp method (http://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method) for details on this algorithm.

The ODE system must be given as the derivative dy/dx = f(x,y) with x in R and y in R^n.

Within this class the function f() is a pure virtual function, which must be reimplemented in a derived class.

Definition at line 45 of file rkodesolver.h.


Constructor & Destructor Documentation

template<typename T, int D>
RkOdeSolver< T, D >::RkOdeSolver ( const T &  x,
const Matrix< T, D, 1 > &  y,
const T &  dx = 0.0,
const T &  eps = 1e-6 
)

Constructor.

Parameters:
xInitial integration parameter
yInitial function values of function to integrate
dxInitial guess for step size. Will be automatically adjusted to guarantee required precision. dx > 0
epsRelative precision. eps > 0.

Initialises the solver with start conditions.

Definition at line 209 of file rkodesolver.h.

   : m_x(x)
{
   Y(y);
   dX(dx);
   Eps(eps);
}

Here is the call graph for this function:

template<typename T , int D>
RkOdeSolver< T, D >::~RkOdeSolver ( void  ) [virtual]

Destructor.

Definition at line 223 of file rkodesolver.h.

{
}

Member Function Documentation

template<typename T , int D>
const T & RkOdeSolver< T, D >::dX ( void  ) const [inline]

Get current estimated step size dX.

Returns:
Reference of dX value.

Definition at line 195 of file rkodesolver.h.

{
   return m_step;
}

Here is the caller graph for this function:

template<typename T, int D>
void RkOdeSolver< T, D >::dX ( const T &  a)

Set estimated step size dX.

Parameters:
aThe value to be set.

Definition at line 231 of file rkodesolver.h.

{
   if (a <= 0.0)
   {
      kError() << "RkOdeSolver: dx must be > 0";
      m_step = 0.001;            // a very arbitrary value
      return;
   }

   m_step = a;
}
template<typename T , int D>
const Matrix< T, D, 1 > & RkOdeSolver< T, D >::dYdX ( void  ) const [inline]

Get current dy/dx value.

Returns:
Reference of dy/dx vector.

Definition at line 188 of file rkodesolver.h.

{
   return m_dydx;
}
template<typename T , int D>
const T & RkOdeSolver< T, D >::Eps ( void  ) const [inline]

Get current presision.

Returns:
Reference of precision value.

Definition at line 202 of file rkodesolver.h.

{
   return m_eps;
}

Here is the caller graph for this function:

template<typename T, int D>
void RkOdeSolver< T, D >::Eps ( const T &  a)

Set estimated presision.

Parameters:
aThe value to be set.

Definition at line 245 of file rkodesolver.h.

{
   if (a <= 0.0)
   {
      kError() << "RkOdeSolver: eps must be > 0";
      m_eps = 1e-5;              // a very arbitrary value
      return;
   }

   m_eps = a;
}
template<typename T, int D>
virtual Matrix<T,D,1> RkOdeSolver< T, D >::f ( const T &  x,
const Matrix< T, D, 1 > &  y 
) const [protected, pure virtual]

ODE function.

Parameters:
xIntegration value
yFunction value
Returns:
Derivation

This purely virtual function returns the value of dy/dx for the given parameter values of x and y. This function is integrated by the Runge-Kutta algorithm.

template<typename T, int D>
void RkOdeSolver< T, D >::integrate ( const T &  dx)

Integrates the ordinary differential equation from the current x value to x+dx.

Parameters:
dxx-interval size to integrate over starting from x. dx may be negative.

The integration is performed by calling rkStepCheck() repeatedly until the desired x value is reached.

Definition at line 261 of file rkodesolver.h.

{
   if (deltaX == 0)
   {
      return;                   // nothing to integrate
   }

   // init dydx
   m_dydx = f(m_x, m_y);

   static const unsigned int maxiter = 10000;
   const T x2 = m_x + deltaX;

   unsigned int iter;
   for (iter=0;
        iter<maxiter && !rkStepCheck(x2-m_x);
        ++iter)
   {
   }

   if (iter > maxiter)
   {
      kWarning()
         << "RkOdeSolver: More than " << maxiter
         << " iterations in RkOdeSolver::integrate" << endl;
   }
}

Here is the caller graph for this function:

template<typename T, int D>
Matrix< T, D, 1 > RkOdeSolver< T, D >::rkStep ( const T &  dx,
Matrix< T, D, 1 > &  yerr 
) const [private]

Perform one Runge-Kutta integration step forward with step size ::m_step.

Parameters:
dxStep size relative to current x value.
yerrReference to vector in which the estimated error made in y is returned.
Returns:
The y value after the step at x+dx.

Stored current x,y values are not adjusted.

Definition at line 388 of file rkodesolver.h.

{
   static const T
      a2=0.2, a3=0.3, a4=0.6, a5=1.0, a6=0.875,
      b21=0.2,
      b31=3.0/40.0,       b32=9.0/40.0,
      b41=0.3,            b42=-0.9,        b43=1.2,
      b51=-11.0/54.0,     b52=2.5,         b53=-70.0/27.0, b54=35.0/27.0,
      b61=1631.0/55296.0, b62=175.0/512.0, b63=575.0/13824.0,
      b64=44275.0/110592.0, b65=253.0/4096.0,
      c1=37.0/378.0, c3=250.0/621.0, c4=125.0/594.0, c6=512.0/1771.0,
      dc1=c1-2825.0/27648.0,  dc3=c3-18575.0/48384.0,
      dc4=c4-13525.0/55296.0, dc5=-277.0/14336.0, dc6=c6-0.25;

   Matrix<T,D,1> ak2 = f(m_x + a2*dx,
                         m_y + dx*b21*m_dydx);             // 2. step
   Matrix<T,D,1> ak3 = f(m_x + a3*dx,
                         m_y + dx*(b31*m_dydx + b32*ak2)); // 3.step
   Matrix<T,D,1> ak4 = f(m_x + a4*dx,
                         m_y + dx*(b41*m_dydx + b42*ak2
                                   + b43*ak3));           // 4.step
   Matrix<T,D,1> ak5 = f(m_x + a5*dx,
                         m_y + dx*(b51*m_dydx + b52*ak2
                                   + b53*ak3 + b54*ak4)); // 5.step
   Matrix<T,D,1> ak6 = f(m_x + a6*dx,
                         m_y + dx*(b61*m_dydx + b62*ak2
                                   + b63*ak3 + b64*ak4
                                   + b65*ak5));           // 6.step
   yerr       = dx*(dc1*m_dydx + dc3*ak3 + dc4*ak4 + dc5*ak5 + dc6*ak6);
   return m_y + dx*( c1*m_dydx +            c3*ak3 +  c4*ak4 +  c6*ak6);
}
template<typename T, int D>
bool RkOdeSolver< T, D >::rkStepCheck ( const T &  dx) [private]

Perform one integration step with a tolerable relative error given by ::mErr.

Parameters:
dxMaximal step size, may be positive or negative depending on integration direction.
Returns:
Flag indicating if made absolute integration step was equal |dx | (true) less than |dx | (false).

A new estimate for the maximum next step size is saved to ::m_step. The new values for x, y and f are saved in ::m_x, ::m_y and ::m_dydx.

Definition at line 294 of file rkodesolver.h.

{
   static const T safety =  0.9;
   static const T pshrnk = -0.25;
   static const T pgrow  = -0.2;

   // reduce step size by no more than a factor 10
   static const T shrinkLimit = 0.1;
   // enlarge step size by no more than a factor 5
   static const T growthLimit = 5.0;
   // errmax_sl = 6561.0
   static const T errmax_sl = std::pow(shrinkLimit/safety, 1.0/pshrnk);
   // errmax_gl = 1.89e-4
   static const T errmax_gl = std::pow(growthLimit/safety, 1.0/pgrow);

   static const unsigned int maxiter = 100;

   if (dx_requested == 0)
   {
      return true;              // integration done
   }

   Matrix<T,D,1> ytmp, yerr, t;

   bool stepSizeWasMaximal;
   T dx;
   if (std::abs(dx_requested) > m_step)
   {
      stepSizeWasMaximal = true;
      dx = dx_requested>0 ? m_step : -m_step;
   }
   else
   {
      stepSizeWasMaximal = false;
      dx = dx_requested;
   }

   // generic scaling factor
   // |y| + |dx * dy/dx| + 1e-15
   Matrix<T,D,1> yscal
      = (m_y.cwise().abs() + (dx*m_dydx).cwise().abs()).cwise()
      + 1e-15;

   unsigned int iter = 0;
   T errmax = 0;
   do
   {
      if (errmax >= 1.0)
      {
         // reduce step size
         dx *= errmax<errmax_sl ? safety * pow(errmax, pshrnk) : shrinkLimit;
         stepSizeWasMaximal = true;
         if (m_x == m_x + dx)
         {
            // stepsize below numerical resolution
            kError()
               << "RkOdeSolver: stepsize underflow in rkStepCheck"
               << endl;
         }
         // new dx -> update scaling vector
         yscal
            = (m_y.cwise().abs()
               + (dx*m_dydx).cwise().abs()).cwise()
            + 1e-15;
      }

      ytmp   = rkStep(dx, yerr); // try to make a step forward
      t      = (yerr.cwise() / yscal).cwise().abs(); // calc the error vector
      errmax = t.maxCoeff()/m_eps;    // calc the rel. maximal error
      ++iter;
   } while ((iter < maxiter) && (errmax >= 1.0));

   if (iter >= maxiter)
   {
      kError()
         << "RkOdeSolver: too many iterations in rkStepCheck()";
   }

   if (stepSizeWasMaximal)
   {
      // estimate next step size if used step size was maximal
      m_step =
         std::abs(dx)
         * (errmax>errmax_gl ? safety * pow(errmax, pgrow) : growthLimit);
   }
   m_x    += dx;                // make step forward
   m_y     = ytmp;              // save new function values
   m_dydx  = f(m_x,m_y);        // and update derivatives

   return (std::abs(dx) < std::abs(dx_requested));
}
template<typename T , int D>
const T & RkOdeSolver< T, D >::X ( void  ) const [inline]

Get current x value.

Returns:
Reference of x value.

Definition at line 160 of file rkodesolver.h.

{
   return m_x;
}
template<typename T, int D>
void RkOdeSolver< T, D >::X ( const T &  a) [inline]

Set current x value.

Parameters:
aThe value to be set.

Definition at line 167 of file rkodesolver.h.

{
   m_x = a;
}
template<typename T , int D>
const Matrix< T, D, 1 > & RkOdeSolver< T, D >::Y ( void  ) const [inline]

Get current y value.

Returns:
Reference of y vector.

Definition at line 174 of file rkodesolver.h.

{
   return m_y;
}

Here is the caller graph for this function:

template<typename T, int D>
void RkOdeSolver< T, D >::Y ( const Matrix< T, D, 1 > &  a) [inline]

Set current y values.

Parameters:
aThe vector to be set.

Definition at line 181 of file rkodesolver.h.

{
   m_y = a;
}

Member Data Documentation

template<typename T, int D>
Matrix<T,D,1> RkOdeSolver< T, D >::m_dydx [private]

current value of dy/dx

Definition at line 148 of file rkodesolver.h.

template<typename T, int D>
T RkOdeSolver< T, D >::m_eps [private]

allowed relative error

Definition at line 151 of file rkodesolver.h.

template<typename T, int D>
T RkOdeSolver< T, D >::m_step [private]

estimated step size for next Runge-Kutta step

Definition at line 153 of file rkodesolver.h.

template<typename T, int D>
T RkOdeSolver< T, D >::m_x [private]

current x value

Definition at line 144 of file rkodesolver.h.

template<typename T, int D>
Matrix<T,D,1> RkOdeSolver< T, D >::m_y [private]

current y value

Definition at line 146 of file rkodesolver.h.


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