numerics/runge_kutta.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_RUNGE_KUTTA_H
#define _MATHTL_RUNGE_KUTTA_H

#include <algebra/vector.h>
#include <algebra/triangular_matrix.h>
#include <numerics/one_step_scheme.h>

namespace MathTL
{
  template <class VECTOR>
  class ExplicitRungeKuttaScheme
    : public OneStepScheme<VECTOR>
  {
  public:
    enum Method {
      RK12,       // s=2, p=2(1) (embedded scheme is Euler scheme)
      RK23,       // s=3, p=3(2) (embedded scheme is the Runge/improved Euler scheme), see [S]
      Fehlberg34, // s=5 (effective: 4), p=4(3), FSAL (uses the classical RK scheme), see [DB]
      DoPri45,    // s=7 (effective: 6), p=5(4), FSAL
      DoPri78     // s=13 (effective: 12), p=8(7), FSAL
    };

    ExplicitRungeKuttaScheme(const Method method);

    void increment(const AbstractIVP<VECTOR>* ivp,
                   const double t_m, const VECTOR& u_m,
                   const double tau,
                   VECTOR& u_mplus1,
                   VECTOR& error_estimate,
                   const double tolerance = 1e-2) const;

    int order() const { return p; }

  protected:
    bool fsal;

    Vector<double> c;

    LowerTriangularMatrix<double> A;

    Vector<double> b, bhat;

    int p;
  };
}

#include "numerics/runge_kutta.cpp"

#endif