A..0.0.15 chapter8/p164.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 164
   2 /* **  EULER, HEUN, AND RUNGE-KUTTA  **
   3    ** METHODS FOR SYSTEM OF ORDINARY **
   4    **    DIFFERENTIAL EQUATIONS-1    **
   5    NFUNC : FUNCTION NUMBER
   6    METHOD : METHOD NUMBER
   7    N : NUMBER OF UNKNOWN FUNCTIONS
   8    T : TIME VARIABLE
   9    X : APPROXIMATE VALUE OF SOLUTION
  10    (T0,X0) : NITIAL VALUE OF (T,X)
  11    H : TIME INCREMENT
  12    TMAX : UPPER LIMIT OF TIME
  13 */
  14 
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 
  18 double func(double t, double x[], double f[], int nfunc);
  19 
  20 int main(void)
  21 {
  22   int nfunc, method, n, i, nstep;
  23   double t0, h, Tmax, t;
  24   double x[10], u[10];
  25   double fk1[10],fk2[10],fk3[10],fk4[10];
  26   printf("LINEAR SYSTEM        : NFUNC=1\n");
  27   printf("VAN DER POL EQUATION :      =2\n");
  28   printf("EXERCISE-2           :      =3\n");
  29   printf("INPUT : NFUNC=\n");
  30   scanf("%d", &nfunc);
  31   printf("EULER       : METHOD=1\n");
  32   printf("HEUN        :       =2\n");
  33   printf("RUNGE-KUTTA :       =3\n");
  34   printf("INPUT : METHOD=\n");
  35   scanf("%d", &method);
  36   printf(" 'INPUT : N=\n");
  37   scanf("%d", &n);
  38   printf(" 'INPUT : T0=\n");
  39   scanf("%lf", &t0);
  40   for (i=0; i<n; i++) {
  41     printf("INPUT : X0(%3d)=\n", i+1);
  42     scanf("%lf", &x[i]);
  43   }
  44   printf(" 'INPUT : H=\n");
  45   scanf("%lf", &h);
  46   printf(" 'INPUT : TMAX=\n");
  47   scanf("%lf", &Tmax);
  48   // CALCULATION
  49   t=t0;
  50   nstep=0;
  51   do {
  52     if (method == 1) {
  53       func(t,x,fk1,nfunc);
  54       for (i=0; i<n; i++)
  55     x[i] += h * fk1[i];
  56     }
  57     else if (method == 2) {
  58       func(t,x,fk1,nfunc);
  59       for (i=0; i<n; i++)
  60     u[i] = x[i] + h * fk1[i];
  61       func(t+h,u,fk2,nfunc);
  62       for (i=0; i<n; i++)
  63     x[i] += h * (fk1[i] + fk2[i]) / 2.0;
  64     }
  65     else {
  66       func(t,x,fk1,nfunc);
  67       for (i=0; i<n; i++)
  68     u[i] = x[i] + h * fk1[i] / 2.0;
  69       func(t+h/2.0,u,fk2,nfunc);
  70       for (i=0; i<n; i++)
  71     u[i] = x[i] + h * fk2[i] / 2.0;
  72       func(t+h/2.0,u,fk3,nfunc);
  73       for (i=0; i<n; i++)
  74     u[i] = x[i] + h * fk3[i];
  75       func(t+h,u,fk4,nfunc);
  76       for (i=0; i<n; i++)
  77     x[i] += h * (fk1[i] + 2.0 * fk2[i] + 2.0 * fk3[i] + fk4[i]) / 6.0;
  78     }
  79     nstep++;
  80     t = t0 + h * nstep;
  81     // OUTPUT
  82     printf("T=%12.4e\n", t);
  83     for (i = 0; i < 5; i++)
  84       printf("  I     X(I)    ");
  85     printf("\n");
  86     for (i=0; i<n; i++) {
  87       printf("%3d%12.4e ", i+1, x[i]);
  88       if ((i+1) % 5 == 0 || (i+1) == n)
  89     printf("\n");
  90     }
  91   }
  92   while ((h > 0.0 && t+h < Tmax+0.1*h) || (h < 0.0 && t+h > Tmax+0.1*h));
  93 }
  94 
  95 double func(double t, double x[], double f[], int nfunc)
  96 {
  97   if (nfunc == 1) {
  98     f[0] = -4.0 * x[0] + 2.0 * x[1];
  99     f[1] = -3.0 * x[0] + x[1];
 100   }
 101   else if (nfunc == 2) {
 102     f[0] = x[1];
 103     f[1] = (1.0 - x[0] * x[0]) * x[1] - x[0];
 104   }
 105   else if (nfunc == 3) {
 106     f[0] = x[1];
 107     f[1] = -x[0] - 0.2 * x[1];
 108     if (x[1] < 0.0) f[1] = f[1] - 1.0;
 109     if (x[1] > 0.0) f[1] = f[1] + 1.0;
 110   }
 111   return 0.0; // ここに来ない
 112 }