A..0.0.5 chapter3/p58.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 58
   2 /* ** FINITE DIFFERENCE METHOD **
   3    **      FOR -D2U/Dx2=F      **
   4    N : NUMBER OF DIVISIONS
   5    H : MESH SIZE
   6    A,B,C : COEFFICIENT MATRIX
   7    F : FREE TERM & SOLUTION VECTOR
   8 */
   9 
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #define NDIM (100)
  13 
  14 void trid(double a[], double b[], double c[], double f[], int n);
  15   
  16 int main(void)
  17 {
  18   int n, i;
  19   double h, h2;
  20   double a[NDIM], b[NDIM], c[NDIM], f[NDIM], x[NDIM];
  21   printf("INPUT : N=\n");
  22   scanf("%d", &n);
  23   if (n <= 1) exit(0);
  24   h = 1.0 / n;
  25   h2 = h * h;
  26   for (i=0; i<n; i++) {
  27     a[i] = -1.0 / h2;
  28     b[i] = 2.0 / h2;
  29     c[i] = -1.0 / h2;
  30     f[i] = 1.0;
  31   }
  32   b[n-1] = 1.0 / h2;
  33   f[n-1] = 0;
  34   trid(a,b,c,f,n);
  35   printf("MESH SIZE=%12.4e\n", h);
  36   printf("** NODAL VALUES OF SOLUTION **\n");
  37   for (i = 0; i < 5; i++)
  38     printf("  I     U(I)    ");
  39   printf("\n");
  40   for (i=0; i<n; i++) {
  41     printf("%3d%12.4e ", i+1, f[i]);
  42     if ((i+1) % 5 == 0 || (i+1) == n)
  43     printf("\n");
  44   }
  45 }
  46 
  47 void trid(double a[], double b[], double c[], double f[], int n)
  48 {
  49   int i;
  50   double aa;
  51   // FORWARD ELIMITNATION
  52   for (i = 0; i < n-1; i++) {
  53     aa = a[i+1] / b[i];
  54     b[i+1] -= aa * c[i];
  55     f[i+1] -= aa * f[i];
  56   }
  57   // BACKWARD SUBSTITUTION
  58   f[n-1] /= b[n-1];
  59   for (i=n-2; i>= 0; i--)
  60     f[i] = (f[i] - c[i] * f[i+1]) / b[i];
  61 }