                                                 {
 //      Take a step using RK4
 
 //      Initialization
         Y0 = Y; Y1 = Y;
 
 //      Step 1
         (*F)(Y1,Yh);                    // slope in the beginning
         Yh *= (0.5*h);   Y1 += Yh;      // Y1 = Y1 + (h/2)*Yh
         Yh *= (1.0/3.0); Y  += Yh;      // Y  = Y  + (h/6)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 2
         (*F)(Y1,Yh);     Y1  = Y0;      // slope in the middle
         Yh *= (0.5*h);   Y1 += Yh;      // Y1 = Y0 + (h/2)*Yh
         Yh *= (2.0/3.0); Y  += Yh;      // Y  = Y  + (h/3)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 3
         (*F)(Y1,Yh);                    // slope in the middle again
         Yh *= h;          Y0 += Yh;     // Y0 = Y0 + h*Yh
         Yh *= (1.0/3.0);  Y  += Yh;     // Y  = Y  + (h/3)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 4
         (*F)(Y0,Yh);                    // slope at the end
         Yh *= (h/6.0);    Y += Yh;      // Y  = Y  + (h/6)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
 
         return Y;
     }

◆ operator()() [2/2]

template<class T>

T & Algorithms::RK4< T >::operator()	(	T &	Y,
		double	h,
		AbstFunctor< T > *	F,
		size_t	dir
	)

Definition at line 667 of file lib-algorithms.h.

                                                             {
 //      Take a step using RK4
 
 //      Initialization
         Y0 = Y; Y1 = Y;
 
 //      Step 1
         (*F)(Y1,Yh,dir);                    // slope in the beginning
         Yh *= (0.5*h);   Y1 += Yh;      // Y1 = Y1 + (h/2)*Yh
         Yh *= (1.0/3.0); Y  += Yh;      // Y  = Y  + (h/6)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 2
         (*F)(Y1,Yh,dir);     Y1  = Y0;      // slope in the middle
         Yh *= (0.5*h);   Y1 += Yh;      // Y1 = Y0 + (h/2)*Yh
         Yh *= (2.0/3.0); Y  += Yh;      // Y  = Y  + (h/3)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 3
         (*F)(Y1,Yh,dir);                    // slope in the middle again
         Yh *= h;          Y0 += Yh;     // Y0 = Y0 + h*Yh
         Yh *= (1.0/3.0);  Y  += Yh;     // Y  = Y  + (h/3)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 //      Step 4
         (*F)(Y0,Yh,dir);                    // slope at the end
         Yh *= (h/6.0);    Y += Yh;      // Y  = Y  + (h/6)*Yh
 //      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
 
         return Y;
     }