ode_runge4.hpp

/* ode_runge4.hpp */
/* $Id$ */
#ifndef __ode_runge4_hpp
#define __ode_runge4_hpp

#include <mth/ode.hpp>

namespace ode
{
  template <typename T, typename X, typename F> class Runge4_1D : public Solver1D<T, X, F>
  {
  public:
    
    bool operator()(const T &t,  
            const T &dt, 
            F &f,        
            X &x,        
            X &dx        
            ) const 
    {return step(t, dt, f, x, dx);}
    
    bool operator()(const T &t,  
            const T &dt, 
            F &f,        
            X &x,        
            X &dx,       
            X &theta     
            ) const
    {return step(t, dt, f, x, dx, theta);}
    
    
    
  protected:
    
    virtual bool step(const T &t,  
              const T &dt, 
              F &f,        
              X &x,        
              X &dx        
              ) const;
    
    bool step(const T &t,  
          const T &dt, 
          F &f,        
          X &x,        
          X &dx,       
          X &theta     
          ) const;
  };
  
  template <typename T, typename X, typename F> class Runge4 : public Solver<T, X, F>
  {
  protected:
    
    virtual bool step(const T &t,              
              const T &dt,             
              const std::vector<F> &f, 
              std::vector<X> &x,       
              std::vector<X> &dx       
              ) const;
  };
  
  
  // ------------------------------------------------------------------------------------------------------------
  // -- Inline implementation -----------------------------------------------------------------------------------
  // ------------------------------------------------------------------------------------------------------------
  
  template <typename T, typename X, typename F>
  bool Runge4_1D<T,X,F>::step(const T &t,  
                  const T &dt, 
                  F &f,        
                  X &x,        
                  X &dx        
                  ) const
  {
    X k1 = f(t, x);
    
    dx = x + k1*dt*(T)0.5; // using dx as temporary variable
    X k2 = f(t + dt*(T)0.5, dx);
    
    dx = x + k2*dt*(T)0.5;
    X k3 = f(t + dt*(T)0.5, dx);
    
    dx = x + k3*dt;
    X k4 = f(t + dt, dx);
    
    dx = (k1 + (T)2*k2 + (T)2*k3 + k4)*dt/(T)6.;
    
    return true;
  }
  
  
  template <typename T, typename X, typename F>
  bool Runge4_1D<T,X,F>::step(const T &t,  
                  const T &dt, 
                  F &f,        
                  X &x,        
                  X &dx,
                  X &theta
                  ) const
  {
    X k1 = f(t, x);
    
    dx = x + k1*dt*(T)0.5; // using dx as temporary variable
    X k2 = f(t + dt*(T)0.5, dx);
    
    dx = x + k2*dt*(T)0.5;
    X k3 = f(t + dt*(T)0.5, dx);
    
    dx = x + k3*dt;
    X k4 = f(t + dt, dx);
    
    theta = (k2 - k3)/(k1 - k2);
    
    dx = (k1 + (T)2*k2 + (T)2*k3 + k4)*dt/(T)6.;
    
    return true;
  }
  
  /* Inline implementation for one equation solver. */
  template <typename T, typename X, typename F>
  bool Runge4<T,X,F>::step(const T &t,             
               const T &dt,             
               const std::vector<F> &f, 
               std::vector<X> &x,       
               std::vector<X> &dx       
               ) const
  {
    // check sizes
    if (f.size() != x.size() || x.size() != dx.size())
      {
    printf ("ode::Runge4::solve: error, different sizes of f(%zu), x(%zu), dx(%zu)\n",
        f.size(), x.size(), dx.size());
      }
       
    // calc current right parts values
    std::vector<X> k1(f.size(), 0.);
    for (size_t i=0; i<f.size(); i++)
      k1[i] = f[i](t, x);

    // shifted vector
    std::vector<X> s(f.size(), 0.);
    for (size_t i=0; i<f.size(); i++)
      s[i] = x[i] + k1[i]*dt*(T)0.5;
    
    // second order koefficients
    std::vector<X> k2(f.size(), 0.);
    for (size_t i=0; i<f.size(); i++)
      k2[i] = f[i](t + dt*(T)0.5, s);
    
    // shifted vector
    for (size_t i=0; i<f.size(); i++)
      s[i] = x[i] + k2[i]*dt*(T)0.5;
    
    // third order koefficients
    std::vector<X> k3(f.size(), 0.);
    for (size_t i=0; i<f.size(); i++)
      k3[i] = f[i](t + dt*(T)0.5, s);
    
    // shifted vector
    for (size_t i=0; i<f.size(); i++)
      s[i] = x[i] + k3[i]*dt;
    
    // fourth order koefficients
    std::vector<X> k4(f.size(), 0.);
    for (size_t i=0; i<f.size(); i++)
      k4[i] = f[i](t + dt, s);
    
    // resulted functions values
    for (size_t i=0; i<f.size(); i++)
      dx[i] = (k1[i] + (T)2*k2[i] + (T)2*k3[i] + k4[i])*dt/(T)6;
   
    return true;
  }
};

#endif /* ode_runge4.hpp */