/* 
************************************************************************
*  runge.c: 常微分方程式の4次Runge-Kutta法による解法
************************************************************************
*/
       
#include <stdio.h>
#include <math.h>

#define MIN 0.0                               /* xの最小 */
#define MAX 1.0                               /* xの最大 */
#define PI 3.14159265358979323846
 
main()
{
   double f(double x, double y);
   
   double x, y, dx, h, dy, yreal, yerror, maxerr;

   double t1, t2, t3,                /* Runge-Kutta法のための */
          k1, k2, k3, k4;            /* 一時変数              */
   
   FILE *output;                     /* runge.datにデータをセーブ  */
   output = fopen("runge.dat","w");

   printf("Input: dx\n");
   scanf("%lf", &dx);
   printf("dx=%g\n", dx);
   
   x = MIN;                                /* 初期条件 */
   y = 0.0;
   yreal = y;
   yerror = 0.0;
   maxerr=0.0;
   
   fprintf(output, "%f\t%f\t%f\t%e\n", x, y, yreal, yerror);

   while(x <= MAX)
   {
      h=dx/2.0;
      k1 = dx*f(x, y);
      t1 = y+0.5*k1;
      k2 = dx*f(x+h, t1);
      t2 = y+0.5*k2;
      k3 = dx*f(x+h, t2);
      t3 = y+k3;
      k4 = dx*f(x + dx, t3);
      x += dx;                              /* xを更新 */    
      y += (k1+2*k2+2*k3+k4)/6.0;           /* yを更新 */       

      yreal=sin(2.0*PI*x);                  /* 解析解 */ 
      yerror=y-yreal;                       /* 誤差=数値解-解析解 */
      if(fabs(yerror) > maxerr){
        maxerr=fabs(yerror);                /* 誤差の絶対値の最大値 */ 
      }
      fprintf(output, "%f\t%f\t%f\t%e\n", x, y, yreal, yerror);
   }
  printf("maxerr =%e\n",maxerr);
   printf("data stored in runge.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

/* dy/dx */
double  f(double x, double y)
{
   return(2.0*PI*cos(2.0*PI*x));               /* dy/dx */
}