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

#define N 2                                   /* 方程式の数 */
#define MIN 0.0                               /* xの最小    */
#define MAX 10.0                              /* yの最大    */
 
main()
{
   void rungen(double x, double y[], double dx);
   double f(double x, double y[], int i);
   
   double x, y[N], dx, h, yreal, yerror, maxerr;
   int j;
   
   FILE *output;                              /* rungen.datにデータをセーブ */
   output = fopen("rungen.dat","w");

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

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

/* Runge-Kutta法 */
void rungen(double x, double y[], double dx)
{
   double f(double x, double y[], int i);
   double h=dx/2.0,
          t1[N], t2[N], t3[N],                /* Runge-Kutta法のための */
          k1[N], k2[N], k3[N],k4[N];          /* 一時変数              */
   int i;
 
   for (i=0; i<N; i++) 
     {
       k1[i] = dx*f(x, y, i);
       t1[i] = y[i]+0.5*k1[i];
     }
   for (i=0; i<N; i++) 
     {
       k2[i] = dx*f(x+h, t1, i);
       t2[i] = y[i]+0.5*k2[i];
     }
   for (i=0; i<N; i++)
     {
       k3[i] = dx*f(x+h, t2, i);
       t3[i] = y[i]+    k3[i];
     }
   for (i=0; i<N; i++) 
     {
       k4[i] = dx*f(x + dx, t3, i);
     }
   for (i=0; i<N; i++) y[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/

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