A..0.0.16 chapter8/p170.c

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