/*
*********************************************************************
*  rungen2.c: 2元連立常微分方程式の4次Runge-Kutta法による解法
************************************************************************
*/

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

#define MIN 0.0                               /* xの最小    */
#define MAX 10.0                              /* yの最大    */
 
main()
{
   double f(double x, double y0, double y1, int i);
   
   double x, y0, y1, dx, h, yreal, yerror, maxerr;
   double t10, t11, t20, t21, t30, t31,       /* Runge-Kutta法のための */
          k10, k11, k20, k21, k30, k31 ,k40, k41;  /* 一時変数        */
   int i;
   int j;
   
   FILE *output;                              /* rungen.datにデータをセーブ */
   output = fopen("rungen.dat","w");

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

   while(x <= MAX)
   {

     k10 = dx*f(x, y0, y1, 0);
     t10 = y0+0.5*k10;
     k11 = dx*f(x, y0, y1, 1);
     t11 = y1+0.5*k11;
     
     k20 = dx*f(x+h, t10, t11, 0);
     t20 = y0+0.5*k20;
     k21 = dx*f(x+h, t10, t11, 1);
     t21 = y1+0.5*k21;
     
     k30 = dx*f(x+h, t20, t21, 0);
     t30 = y0+    k30;
     k31 = dx*f(x+h, t20, t21, 1);
     t31 = y1+    k31;
     
     k40 = dx*f(x + dx, t30, t31, 0);
     k41 = dx*f(x + dx, t30, t31, 1);
     
     y0 += (k10+2*k20+2*k30+k40)/6.0;
     y1 += (k11+2*k21+2*k31+k41)/6.0;

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

/* definition of equations - this is the harmonic oscillator */
double  f(double x, double y0, double y1, int i)
{
   if (i == 0) return(y1);               /* d(y0)/dx */
   if (i == 1) return(-y0);              /* d(y1)/dx */
}