/* 
************************************************************************
*  euler.c: Euler法による常微分方程式の解法
************************************************************************
*/

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

#define MIN 0.0                               /* xの最小 */
#define MAX 1.0                               /* xの最大 */
#define PI 3.14159265358979323846
/* #define PI M_PI */

main()
{
  double f(double x, double y);
   
  double x, y, dx, dy, yreal, yerror, maxerr;
   
  FILE *output;                           /* euler.datにデータをセーブ */
  output = fopen("euler.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); */
  fprintf(output, "%e\t%e\t%e\t%e\n", x, y, yreal, yerror);
  
  
  while(x <= MAX)                          
    {
      dy = dx*f(x, y);                      /* yの増分 */  
      x += dx;
      y += dy;                              /* 数値解 */ 
      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); */
      fprintf(output, "%e\t%e\t%e\t%e\n", x, y, yreal, yerror);
      }
  printf("maxerr =%e\n",maxerr);
  printf("data stored in euler.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 */
}