A..0.0.6 chapter3/p65.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 65
   2 /* ** FINITE DIFFERENCE METHOD **
   3    ** FOR -FNU D2U/Dx2+DU/DX=F **
   4    FNU : DIFFUSIVITY
   5    N : NO. OF DIVISIONS
   6    METHOD : METHOD NUMBER
   7    H : MESH SIZE
   8    A,B,C : COEFFICIENT MATRIX
   9    F : FREE TERM & SOLUTION VECTOR
  10 */
  11 
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #define NDIM (100)
  15 void trid(double a[], double b[], double c[], double f[], int n);
  16   
  17 int main(void)
  18 {
  19   int n, i, method;
  20   double fnu, h, h2; // fnu は nu でも良いかも。
  21   double a[NDIM], b[NDIM], c[NDIM], f[NDIM], x[NDIM];
  22   printf("INPUT : FNU=\n");
  23   scanf("%lf", &fnu);
  24   printf("INPUT : N=\n");
  25   scanf("%d", &n);
  26   if (n <= 1) exit(1);
  27   printf("BACKWARD : METHOD=1, CENTRAL : =2, FORWARD : =3\n");
  28   printf("INPUT : METHOD=\n");
  29   scanf("%d", &method);
  30   h = 1.0 / n; // (double) n とする必要もないでしょう。
  31   h2 = h * h;  // Cに冪乗演算子はないので、乗算で代用
  32   for (i=0; i<n-1; i++) {
  33     a[i] = - fnu / h2;
  34     b[i] = 2.0 * fnu / h2;
  35     c[i] = - fnu / h2;
  36     if (method == 1) {
  37       // BACKWARD
  38       a[i] -= 1.0 / h;
  39       b[i] += 1.0 / h;
  40     }
  41     else if (method == 2) {
  42       // CENTRAL
  43       a[i] -= 0.5 / h;
  44       c[i] += 0.5 / h;
  45     }
  46     else {
  47       // FORWARD
  48       b[i] -= 1.0 / h;
  49       c[i] += 1.0 / h;
  50     }
  51     f[i] = 1.0;
  52   }
  53   trid(a,b,c,f,n-1);
  54   printf("DIFFUSIVITY=%12.4e   MESH SIZE=%12.4e   METHOD=%2d\n",
  55      fnu, h, method);
  56   printf("** NODAL VALUES OF SOLUTION **\n");
  57   for (i = 0; i < 5; i++)
  58     printf("  I     U(I)    ");
  59   printf("\n");
  60   for (i=0; i<n-1; i++) {
  61     printf("%3d%12.4e ", i+1, f[i]);
  62     if ((i+1) % 5 == 0 || (i+1) == n-1)
  63     printf("\n");
  64   }
  65 }
  66 
  67 void trid(double a[], double b[], double c[], double f[], int n)
  68 {
  69   int i;
  70   double aa;
  71   // FORWARD ELIMITNATION
  72   for (i=0; i<n-1; i++) {
  73     aa = a[i+1] / b[i];
  74     b[i+1] -= aa * c[i];
  75     f[i+1] -= aa * f[i];
  76   }
  77   // BACKWARD SUBSTITUTION
  78   f[n-1] /= b[n-1];
  79   for (i=n-2; i>= 0; i--)
  80     f[i] = (f[i] - c[i]*f[i+1]) / b[i];
  81 }