// rk4th_fourth.c   fourth order ODE  d^4U(x)/d^4x = F4(x)
//             let y = sin(x)  the initial conditions y, yp, ypp, yppp at x=0
// y = sin(0) = 0.0,  y'=yp=sin'(0)=1.0  y''=ypp=sin''(0) = 0.0
// y'''=yppp=sin'''(0)=-1.0    thus F4(x) = sin(x) for 4th derivative
//

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

double F4(double x, double y)  // y may be unused
{ // return information on 4rd derivative
  return sin(x);
}

int main(int argc, char * argv[])
{
  double k1, k2, k3, k4; // for y, based on yp
  double dk1, dk2, dk3, dk4; // for yp, based on ypp
  double ddk1, ddk2, ddk3, ddk4; // for ypp, based on yppp
  double dddk1, dddk2, dddk3, dddk4; // for yppp, based on F4 RHS
  double x, h, y, yp, ypp, yppp;
  int i;

   printf("Runge-Kutta 4th, on fourth order equation \n");
   x = 0.0;
   y = 0.0;
   yp = 1.0;
   ypp = 0.0;
   yppp = -1.0;
   h = M_PI/32.0;
   printf("given y''''=sin(x) solution is y(x)=sin(x) \n");
   printf("\n");
   printf("             y= sin(x)    y'= cos(x)    y''=-sin(x)    y'''=-cos(x)\n");
   printf("x=%9.6f, y=%9.6f, yp=%9.6f, ypp=%9.6f, yppp=%9.6f, h=%9.6g \n",
          x, y, yp, ypp, yppp, h);
   for(i=0;i<64;i++)
   {
     dddk1 = h*F4(x, y);
     dddk2 = h*F4(x+h/2.0, y+dddk1/2.0);
     dddk3 = h*F4(x+h/2.0, y+dddk2/2.0);
     dddk4 = h*F4(x+h, y+dddk3 );

     ddk1 = h*yppp;
     ddk2 = h*(yppp + dddk1/2.0);
     ddk3 = h*(yppp + dddk2/2.0);
     ddk4 = h*(yppp + dddk3);

     dk1 = h*ypp;
     dk2 = h*(ypp + ddk1/2.0);
     dk3 = h*(ypp + ddk2/2.0);
     dk4 = h*(ypp + ddk3);

     k1 = h*yp;
     k2 = h*(yp + dk1/2.0);
     k3 = h*(yp + dk2/2.0);
     k4 = h*(yp + dk3);

     yppp = yppp + (dddk1 + 2.0*dddk2 + 2.0*dddk3 + dddk4)/6.0;
     ypp =  ypp  + (ddk1  + 2.0*ddk2  + 2.0*ddk3  + ddk4) /6.0;
     yp =   yp   + (dk1   + 2.0*dk2   + 2.0*dk3   + dk4)  /6.0;
     y =    y    + (k1    + 2.0*k2    + 2.0*k3    + k4)   /6.0;
     x = x + h;
     printf("x=%9.6f, y=%9.6f, yp=%9.6f, ypp=%9.6f, yppp=%9.6f \n",
               x, y, yp, ypp, yppp);
   }
   printf("y=sin(x) err    = %g \n", y-sin(x));
   printf("yp=cos(x) err   = %g \n", yp-cos(x));
   printf("ypp=-sin(x) err = %g \n", ypp+sin(x));
   printf("yppp=-cos(x) err= %g \n", yppp+cos(x));

   printf("\n");
   printf("rk4th_fourth.c finished \n");
   return 0;
}