include/apfel/ode.h

0001 //
0002 // APFEL++ 2017
0003 //
0004 // Author: Valerio Bertone: valerio.bertone@cern.ch
0005 //
0006
0007 #pragma once
0008
0009 namespace apfel
0010 {
0011   /**
0012    * @name Runge-Kutta (RK) ODE solvers.
0013    * These functions solve the ordinary differential equation (ODE):
0014    *
0015    *    dy / dt = f(t,y)
0016    *
0017    * where:
0018    *
0019    *    dy = rk4(f(t,y))
0020    *
0021    * so differentiation between lower and upper:
0022    *
0023    *    y += dy(t,y,dt)
0024    *
0025    * U is the type of the 'y' object.
0026    */
0027   ///@{
0028   /**
0029    * @brief Template function that implements the fourth order RK
0030    * algorithm.
0031    * @param f: the function on the r.h.s. of the ODE
0032    * @return the function tha returns the step
0033    */
0034   template<class U>
0035   std::function<U(double const&, U const&, double const&)>
0036   rk4(std::function<U(double const& t, U const& Obj)> const& f)
0037   {
0038     // *INDENT-OFF*
0039     return
0040       [       f            ] (double const& t, U const& y,  double const& dt) -> U{ return
0041       [t,y,dt,f            ] (                 U const& dy1                 ) -> U{ return
0042       [t,y,dt,f,dy1        ] (                 U const& dy2                 ) -> U{ return
0043       [t,y,dt,f,dy1,dy2    ] (                 U const& dy3                 ) -> U{ return
0044       [       f,dy1,dy2,dy3] (                 U const& dy4                 ) -> U{ return
0045       ( dy1 + 2 * dy2 + 2 * dy3 + dy4 ) / 6  ;} (
0046       dt * f( t + dt    , y + dy3     )     );} (
0047       dt * f( t + dt / 2, y + dy2 / 2 )     );} (
0048       dt * f( t + dt / 2, y + dy1 / 2 )     );} (
0049       dt * f( t         , y           )     );} ;
0050     // *INDENT-ON*
0051   }
0052
0053   /**
0054    * @brief Template function that implements the first order RK
0055    * algorithm.
0056    * @param f: the function on the r.h.s. of the ODE
0057    * @return the function tha returns the step
0058    */
0059   template<class U>
0060   std::function<U(double const&, U const&, double const&)>
0061   rk1(std::function<U(double const& t, U const& Obj)> const& f)
0062   {
0063     return [f] (double const& t, U const& y,  double const& dt) -> U{ return dt * f(t, y); } ;
0064   }
0065   ///@}
0066 }