A..0.0.8 chapter3/p75.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 75
   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 : FREETERM & SOLUTION VECTOR
  10 */
  11 
  12 #include <stdio.h>
  13 #include <stdlib.h> // exit()
  14 #define NDIM (100)
  15 
  16 void trid(double a[], double b[], double c[], double f[], int n);
  17   
  18 int main(void)
  19 {
  20   int n, method, i;
  21   double fnu, h, h2;
  22   double a[NDIM], b[NDIM], c[NDIM], f[NDIM];
  23   printf("INPUT : FNU=\n");
  24   scanf("%lf", &fnu);
  25   printf("INPUT : N=\n");
  26   scanf("%d", &n);
  27   if (n <= 1) exit(1);
  28   printf("BACKWARD : METHOD=1,CENTRAL : =2,FORWARD : =3\n");
  29   printf("INPUT : METHOD=\n");
  30   scanf("%d", &method);
  31   h = 1.0 / n;
  32   h2 = h * h;
  33   for (i = 0; i < n-1; i++) {
  34     a[i] = - fnu / h2;
  35     b[i] = 2.0 * fnu / h2;
  36     c[i] = - fnu / h2;
  37     if (method == 1) {
  38       // BACKWARD
  39       a[i] -= 1.0 / h;
  40       b[i] += 1.0 / h;
  41     }
  42     else if (method == 2) {
  43       // CENTRAL
  44       a[i] -= 0.5 / h;
  45       c[i] += 0.5 / h;
  46     }
  47     else {
  48       // FORWARD
  49       b[i] -= 1.0 / h;
  50       c[i] += 1.0 / h;
  51     }
  52     f[i] = 1.0;
  53   }
  54   a[n-1] = - fnu / h2;
  55   b[n-1] = fnu / h2;
  56   f[n-1] = 0.0;
  57   if (method == 2) f[n-1] = 0.5;
  58   if (method == 3) f[n-1] = 1.0;
  59   trid(a, b, c, f, n);
  60   printf("DIFFUSIVITY=%12.4e   MESH SIZE=%12.4e   METHOD=%2d\n",
  61      fnu, h, method);
  62   printf("** NODAL VALUES OF SOLUTION **\n");
  63   for (i = 0; i < 5; i++)
  64     printf("  I     U(I)    ");
  65   printf("\n");
  66   for (i=0; i<n; i++) {
  67     printf("%3d%12.4e ", i+1, f[i]);
  68     if ((i+1) % 5 == 0 || (i+1) == n)
  69       printf("\n");
  70   }
  71 }
  72 
  73 void trid(double a[], double b[], double c[], double f[], int n)
  74 {
  75   int i;
  76   double aa;
  77   // FORWARD ELIMINATION
  78   for (i=0; i < n-1; i++) {
  79     aa = a[i+1] / b[i];
  80     b[i+1] -= aa * c[i];
  81     f[i+1] -= aa * f[i];
  82   }
  83   // BACKWARD SUBSTITUTION
  84   f[n-1] /= b[n-1];
  85   for (i = n-2; i >= 0; i--)
  86     f[i] = (f[i]-c[i]*f[i+1]) / b[i];
  87 }