gistfile1.txt

I would like to test numerical recipe's version 3 (nm3) implementation of runge kutta 4 (rk4).

The example ODE that I would like to test is: dy/dx = 3yx^(2)y; initial conditions: x_0=1, y_0=2
This is a separable ODE and therefore the solution is y = (2/e)*e^(x^3)  or 2*e^[(x^3)-1] 
 
I would like to pick h = 0.1
 
   Therefore, via rk4:

y_1 = y_0 + h* Fourth-Order-Taylor(x_0, y_0, h)
y_1 =  2  + (0.1/6)(k_1 + 2k_2 + 2k_3 + k_4)
 
   Calculating k_n:
 
            k_1 = f[ x_0 , y_0 ] = 3(1^2)(2) = 6
            k_2 = f[ (x_0 + h/2) , (y_0 + h/2)*k_1 ] = f[1.05 , 2] = 7.607
            k_3 = ....... = 2.872
            k_4 = ....... = 10.117
   
   Calculating y_1:

            y_1 = 2 + (0.1/6) * ( 6 + 2(7.607) + 2(2.872) + 10.117 ) 
            y_1 = 2.787    and x_1 = x_0 + h = 1+ 0.1 = 1.1 

------------------------------------------------------------------------------------------------

#include <iostream>
#include <nr3.h>
#include <rk4.h>

using namespace std;

void derivs(const Doub x, VecDoub_I & y, VecDoub_O & dydx)
{
    dydx[0]=6;
    dydx[1]=7.607
}

int main ()
{
	VecDoub y(1),dydx(1);
	Doub x, xmin, xmax, kmax = 300000, h = 0.0001;
	VecDoub yout(1);
	int k;

	xmin=0; xmax=10000;

	y[0]=2;
	y[1]=2

	derivs(xmin,y,dydx);

	for(k=0; k < kmax; k++)
	{
	   x = xmin + k*h;

	   rk4(y, dydx, x, h, yout, derivs);
	   
           cout << x << " " << yout[0] << " " << yout[1] << endl;
	   
           y[0]=yout[0];
	   y[1]=yout[1];
	   
           derivs(x,y,dydx);
        }

return 0;

}