/* rk4th_second.c   demo fourth order Runge-Kutta on second order */
/*  given  d^2y/dx^2 = F(x, y, y')    solve for y(x)              */
/*  above form can be from A y'' + B y' + C y = G(x, y, y')       */
/*  where  F(x, y, y') = (G(x, y, y') - B y' - C y)/A             */

#include <stdio.h>
#include <math.h>

#define F(x,y,yp)  (-sin(x))
#define G(x)       (-x*x*x*x+6.0*x*x*x+21.0*x*x+20.0*x+9.0)
#define FG(x,y,yp) (G(x)-2.0*yp+y)

int main(int argc, char * argv[])
{
  float k1, k2, k3, k4; /* for first derivative yp, based on F */
  float dk1, dk2, dk3, dk4; /* for solution y, based on yp */
  float x, h, y, yp;
  int i;

   printf("\n Runge-Kutta 4th, second order equation \n");
   x = 0.0;
   y = 0.0;
   yp = 1.0;
   h = M_PI/32.0;
   printf("given y''=-sin(x) solution is y(x)=sin(x) \n");
   printf("using second order RK 4th \n");
   printf("x=%9.6f, y=%9.6f, yp=%9.6f, err=%g \n",
          x, y, yp, y-sin(x));
   for(i=0;i<32;i++)
   {
     k1 =  h*F(x, y, yp);
     dk1 = h*yp;
     k2 =  h*F(x+h/2.0, y+dk1/2.0, yp+k1/2.0);
     dk2 = h*(yp+k1/2.0);
     k3 =  h*F(x+h/2.0, y+dk2/2.0, yp+k2/2.0);
     dk3 = h*(yp+k2/2.0);
     k4 =  h*F(x+h, y+dk3, yp+k3);
     dk4 = h*(yp+k3);
     yp =  yp + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
     y  =  y + (dk1 + 2.0*dk2 + 2.0*dk3 + dk4)/6.0;
     x =   x + h;
     printf("x=%9.6f, y=%9.6f, yp=%9.6f, err=%g \n",
            x, y, yp, y-sin(x));
   }

   printf("\n");
   x = 0.0;
   y = 0.0;
   yp = 1.0;
   h = M_PI/32.0;
   printf("given y''=-sin(x) solution is y(x)=sin(x) \n");
   printf("using only first order RK 4th \n");
   printf("x=%9.6f, y=%9.6f, yp=%9.6f, err=%g \n",
          x, y, yp, y-sin(x));
   for(i=0;i<32;i++)
   {
     k1 =  h*F(x, y, yp);
     k2 =  h*F(x+h/2.0, y+h*k1/2.0, yp+k1/2.0);
     k3 =  h*F(x+h/2.0, y+h*k2/2.0, yp+k2/2.0);
     k4 =  h*F(x+h, y+h*k3, yp+k3);
     yp =  yp + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
     y  =  y + h*yp;
     x =   x + h;
     printf("x=%9.6f, y=%9.6f, yp=%9.6f, err=%g \n",
            x, y, yp, y-sin(x));
   }

   printf("\n");
   x = 0.0;
   y = 5.0;
   yp = 4.0;
   h = 0.01;
   printf("given y'' +2y' +y = G(x) \n");
   printf("G(x)=-x^4+6x^3+21x^2+20x+9 \n");
   printf("Solution y(x)=x^4+2x^3+3x^2+4x+5 \n");
   printf("x=%9.6f, y=%9.6f, yp=%9.6f, yerr=%g, yperr=%g \n",
          x, y, yp, y-(x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0),
          yp-(4.0*x*x*x+6.0*x*x+6.0*x+4.0));
   for(i=0;i<100;i++)
   {
     k1 =  h*F(x, y, yp);
     dk1 = h*yp;
     k2 =  h*F(x+h/2.0, y+dk1/2.0, yp+k1/2.0);
     dk2 = h*(yp+k1/2.0);
     k3 =  h*F(x+h/2.0, y+dk2/2.0, yp+k2/2.0);
     dk3 = h*(yp+k2/2.0);
     k4 =  h*F(x+h, y+dk3, yp+k3);
     dk4 = h*(yp+k3);
     yp =  yp + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
     y  =  y + (dk1 + 2.0*dk2 + 2.0*dk3 + dk4)/6.0;
     x =   x + h;
     printf("x=%9.6f, y=%9.6f, yp=%9.6f, yerr=%g, yperr=%g \n",
            x, y, yp, y-(x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0),
            yp-(4.0*x*x*x+6.0*x*x+6.0*x+4.0));
   }

   printf("\n");
   x = 0.0;
   y = 5.0;
   yp = 4.0;
   h = 0.01;
   printf("given y'' +2y' +y = G(x) \n");
   printf("G(x)=-x^4+6x^3+21x^2+20x+9 \n");
   printf("Solution y(x)=x^4+2x^3+3x^2+4x+5 \n");
   printf("x=%9.6f, y=%9.6f, yp=%9.6f, yerr=%g, yperr=%g \n",
          x, y, yp, y-(x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0),
          yp-(4.0*x*x*x+6.0*x*x+6.0*x+4.0));
   for(i=0;i<100;i++)
   {
     k1 =  h*F(x, y, yp);
     k2 =  h*F(x+h/2.0, y+h*k1/2.0, yp+k1/2.0);
     k3 =  h*F(x+h/2.0, y+h*k2/2.0, yp+k2/2.0);
     k4 =  h*F(x+h, y+h*k3, yp+k3);
     yp =  yp + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
     y  =  y + h*yp;
     x =   x + h;
     printf("x=%9.6f, y=%9.6f, yp=%9.6f, yerr=%g, yperr=%g \n",
            x, y, yp, y-(x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0),
            yp-(4.0*x*x*x+6.0*x*x+6.0*x+4.0));
   }
   return 0;
}