A..0.0.3 chapter1/p29.c

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