A..0.0.2 chapter1/p20.c

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