runge_kutta.hpp

View page source

C++ runge_kutta Source Code

# include <cmpad/vector.hpp>

namespace cmpad { // BEGIN cmpad namespace

// BEGIN PROTOTYPE
template <class Vector, class Fun>
Vector runge_kutta(
   const Fun&                           fun ,
   const Vector&                        yi  ,
   const typename Vector::value_type    tf  ,
   size_t                               ns  )
// END PROTOTYPE
{
   // scalar_type
   typedef typename Vector::value_type scalar_type;
   //
   // two, six
   scalar_type two = scalar_type(2.0);
   scalar_type six = scalar_type(6.0);
   // n
   size_t n = yi.size();
   //
   // h
   scalar_type h  = tf / scalar_type( double(ns) );
   //
   // k1, k2, k3, k4, y_tmp
   Vector k1(n), k2(n), k3(n), k4(n), y_tmp(n);
   //
   //
   // i_step, yf
   Vector yf = yi;
   for(size_t i_step = 0; i_step < ns; ++i_step)
   {  //
      // k1
      k1 = fun(yf);
      //
      // k2
      for(size_t i = 0; i < n; ++i)
         y_tmp[i] = yf[i] + h * k1[i] / two;
      k2 = fun(y_tmp);
      //
      // k3
      for(size_t i = 0; i < n; ++i)
         y_tmp[i] = yf[i] + h * k2[i] / two;
      k3 = fun(y_tmp);
      //
      // k4
      for(size_t i = 0; i < n; ++i)
         y_tmp[i] = yf[i] + h * k3[i];
      k4 = fun(y_tmp);
      //
      // yf
      for(size_t i = 0; i < n; ++i)
         yf[i] = yf[i] + h * (k1[i] + two * k2[i] + two * k3[i] + k4[i]) / six;
   }
   return yf;
}

} // END cmpad namespace