A..0.0.14 chapter7/p146.c

   1 // C  菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 146
   2 /* ** EULER,HEUN,AND RUNGE-KUTTA METHODS **
   3    ** FOR ORDINARY DIFFERENTIAL EQUATIONS **
   4    NFUNC : FUNCTION NUMBER
   5    METHOD : METHOD NUMBER
   6    T : TIME VARIABLE
   7    X : APPROXIMATE VALUE OF SOLUTION
   8    (T0,X0) : INITIAL VALUE OF (T,X)
   9    H : TIME INCREMENT
  10    TMAX : UPPER LIMIT OF TIME
  11 */
  12 
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 
  17 double f(double t, double x, int nfunc);
  18 
  19 int main(void)
  20 {
  21   int nfunc, method, nstep;
  22   double t0, x0, h, Tmax, t, x, fk1, fk2, fk3, fk4;
  23   printf("INPUT : NFUNC=\n");
  24   scanf("%d", &nfunc);
  25   printf("EULER : METHOD=1\n");
  26   printf("HEUN  :       =2\n");
  27   printf("RUNGE-KUTTA   =3\n");
  28   printf("INPUT : METHOD=");
  29   scanf("%d", &method);
  30   printf("INPUT : T0=\n");
  31   scanf("%lf", &t0);
  32   printf("INPUT : x0=\n");
  33   scanf("%lf", &x0);
  34   printf("INPUT : h=\n");
  35   scanf("%lf", &h);
  36   printf("INPUT : Tmax=\n");
  37   scanf("%lf", &Tmax);
  38   t=t0;
  39   x=x0;
  40   nstep=0;
  41   do {
  42     if (method == 1) {
  43       x += h*f(t,x,nfunc);
  44     }
  45     else if (method == 2) {
  46       fk1=f(t,x,nfunc);
  47       fk2=f(t+h,x+h*fk1,nfunc);
  48       x += h*(fk1+fk2)/2.0;
  49     }
  50     else {
  51       fk1=f(t,x,nfunc);
  52       fk2=f(t+h/2.0,x+h*fk1/2.0,nfunc);
  53       fk3=f(t+h/2.0,x+h*fk2/2.0,nfunc);
  54       fk4=f(t+h,x+h*fk3,nfunc);
  55       x += h*(fk1+2.0*fk2+2.0*fk3+fk4)/6.0;
  56     }
  57     nstep++;
  58     t=t0+h*nstep;
  59     printf("T=%12.4e X=%12.4e\n", t, x);
  60   }
  61   while ((h > 0.0 && t+h < Tmax+0.1*h) || (h < 0.0 && t+h > Tmax+0.1*h));
  62   return 0;
  63 }
  64 
  65 double f(double t, double x, int nfunc)
  66 {
  67   if (nfunc == 1) return x;
  68   if (nfunc == 2) return 1.0 - x * x;
  69   if (nfunc == 3) return x * x;
  70   return 0.0;
  71 }