#ifndef _utl_RK4ODEIntegrator_h_ #define _utl_RK4ODEIntegrator_h_ namespace utl { /** \class RK4Iterator RK4ODEIntegrator.h "utl/RK4ODEIntegrator.h" This iterator class will advance the solution for the current point. See class RK4ODEIntegrator for details. \author Darko Veberic \date 20 Jun 2008 \version $Id: RK4ODEIntegrator.h 14717 2009-09-17 20:24:36Z lukas $ \ingroup math */ template class RK4Iterator { public: RK4Iterator(const double x, const Vector& y, DerivativeFunctor& d) : fX(x), fDerivativeFunctor(d) { Copy(d, y, fY); } double GetX() const { return fX; } Vector& GetY() { return fY; } const Vector& GetY() const { return fY; } void Advance(const double dx, const Vector& dYdX) { const unsigned int n = fDerivativeFunctor; // first step: x -> x + dx/2 const double dx2 = 0.5 * dx; Vector y1; for (unsigned int i = 0; i < n; ++i) y1[i] = fY[i] + dx2 * dYdX[i]; // second step const double x1 = fX + dx2; Vector dYdX1; fDerivativeFunctor(x1, y1, dYdX1); for (unsigned int i = 0; i < n; ++i) y1[i] = fY[i] + dx2 * dYdX1[i]; // third step Vector dYdX2; fDerivativeFunctor(x1, y1, dYdX2); for (unsigned int i = 0; i < n; ++i) { y1[i] = fY[i] + dx * dYdX2[i]; dYdX2[i] += dYdX1[i]; } // fourth step fX += dx; fDerivativeFunctor(fX, y1, dYdX1); const double dx6 = (1./6) * dx; for (unsigned int i = 0; i < n; ++i) fY[i] += dx6 * (dYdX[i] + dYdX1[i] + 2.*dYdX2[i]); } void Advance(const double dx) { Vector dYdX; fDerivativeFunctor(fX, fY, dYdX); Advance(dx, dYdX); } RK4Iterator& operator+=(const double dx) { Advance(dx); return *this; } RK4Iterator& operator-=(const double dx) { Advance(-dx); return *this; } private: static void Copy(const unsigned int n, const Vector& a, Vector& b) { for (unsigned int i = 0; i < n; ++i) b[i] = a[i]; } double fX; Vector fY; DerivativeFunctor& fDerivativeFunctor; }; /** \class RK4ODEIntegrator RK4ODEIntegrator.h "utl/RK4ODEIntegrator.h" This class prepares the integration of the Ordinary Differential Equations (ODE) of the form \f[ \frac{{\rm d}y_i}{{\rm d}x} = f_i(x, y_0, \ldots, y_{n-1}) \f] for \f$i\in{0, \ldots, n-1}\f$. ODE is defined through an function class with operator() defined (usualy such classes are called functor in C++). This functor has to define the operator() and unsigned int cast operator: \code class MyODE { public: void operator()(const double x, const double y[10], double dydx[10]) { dydx[0] = ... dydx[1] = ... } operator unsigned int() { return 10; } }; \endcode Operator() calculates derivatives \f${\rm d}y_i/{\rm d}x\f$ of variables \f$y_i\f$ with respect to the independent variable x. The cast operator unsigned int() returns number of OD equations (in this example 10). The class RK4ODEIntegrator returns an iterator that can be advanced with operator += dx for some step dx in independent variable. \code typedef double Vector[10]; RK4ODEIntegrator ode(MyODE()); RK4Iterator it = ode.Begin(x0, y0); it += dx; \endcode Note that the suitable candidates for template parameter Vector are vector classes with fixed size (std::vector<> is not appropriate since the default constructor leaves the object with zero size) that support the operator[] for element access. You can use the C-style arrays double[n] or Offline class utl::SVector. See also utl::RK5ODEIntegrator and testODEIntegrator.cc for a real examples. \author Darko Veberic \date 20 Jun 2008 \version $Id: RK4ODEIntegrator.h 14717 2009-09-17 20:24:36Z lukas $ \ingroup math */ template class RK4ODEIntegrator { public: RK4ODEIntegrator(DerivativeFunctor& d) : fDerivativeFunctor(d) { } template RK4Iterator Begin(const double x, const Vector& y) { return RK4Iterator(x, y, fDerivativeFunctor); } DerivativeFunctor& GetDerivativeFunctor() { return fDerivativeFunctor; } private: DerivativeFunctor& fDerivativeFunctor; }; } #endif