A..0.0.4 chapter2/p39.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 39
   2 /* ** JACOBI OR SOR METHOD **
   3    N : NUMBER OF UNKNOWS
   4    A : LOWER DIAGONAL
   5    F : KNOWN VECTOR & SOLUTION VECTOR
   6    X : SOLUTION VECTOR
   7    N : MAXIMUM NO. OF ITERATIONS
   8    EPS : THRESHOLD VALUE FOR ERROR
   9    METHOD : METHOD NUMBER
  10    OMEGA : RELAXATION PARAMETER
  11 */
  12 
  13 #include <stdio.h>
  14 #include <stdlib.h> // exit()
  15 #include <math.h>   // fabs()
  16 #define NDIM (100)
  17 
  18 void iter(int ndim,
  19       double a[][ndim], double f[], double x[], double xnew[],
  20       double eps, double omega, int n, int m, int method, int *niter);
  21 
  22 int main(void)
  23 {
  24   int n, m, method, i,j,niter;
  25   double eps, omega;
  26   double a[NDIM][NDIM], f[NDIM], x[NDIM], xnew[NDIM];
  27   printf("INPUT : N=\n");
  28   scanf("%d", &n);
  29   if (n <= 1) exit(1);
  30   printf("INPUT : M=\n");
  31   scanf("%d", &m);
  32   printf("INPUT : EPS=\n");
  33   scanf("%lf", &eps);
  34   printf("JACOBI : METHOD=0, SOR : =1\n");
  35   printf("INPUT : METHOD=\n");
  36   scanf("%d", &method);
  37   omega=1.0;
  38   if (method == 1) {
  39     printf("INPUT : OMEGA=\n");
  40     scanf("%lf", &omega);
  41   }
  42   for (i=0; i<n; i++)
  43     for (j=0; j<n; j++) {
  44       printf("INPUT : A(%3d,%3d)=\n", i+1, j+1);
  45       scanf("%lf", &a[i][j]);
  46     }
  47   for (i=0; i<n; i++) {
  48     printf("INPUT : F(%3d)=\n", i+1);
  49     scanf("%lf", &f[i]);
  50   }
  51   for (i=0; i<n; i++) {
  52     printf("INPUT : X(%3d)=\n", i+1);
  53     scanf("%lf", &x[i]);
  54   }
  55   // 試しに表示する
  56   printf("n=%d, m=%d, eps=%f, method=%d\n", n, m, eps, method);
  57   for (i=0; i<n; i++) {
  58     for (j=0; j<n; j++)
  59       printf("%f ", a[i][j]);
  60     printf("   %f %f", f[i], x[i]);
  61     printf("\n");
  62   }
  63   iter(NDIM,a,f,x,xnew,eps,omega,n,m,method,&niter);
  64   printf("NO. OF ITERATIONS=%3d   EPS=%12.4e\n", niter, eps);
  65   for (i = 0; i < 5; i++)
  66     printf("  I     X(I)    ");
  67   printf("\n");
  68   for (i=0; i<n; i++) {
  69     printf("%3d%12.4e ", i+1, x[i]);
  70     if ((i+1) % 5 == 0 || (i+1) == n)
  71       printf("\n");
  72   }
  73 }
  74 
  75 void iter(int ndim,
  76       double a[][ndim], double f[], double x[], double xnew[],
  77       double eps, double omega, int n, int m, int method, int *niter)
  78 {
  79   int i, j;
  80   double error, xmax, xtemp, dx;
  81   *niter = 0;
  82   do {
  83     error = 0.0;
  84     xmax = 0.0;
  85     for (i=0; i<n; i++) {
  86       xtemp = f[i];
  87       for (j=0; j<n; j++)
  88     if (j != i) xtemp -= a[i][j] * x[j];
  89       xtemp /= a[i][i];
  90       if (method == 1)
  91     xtemp = omega * xtemp + (1.0 - omega) * x[i];
  92       dx = fabs(xtemp - x[i]);
  93       if (dx > error) error = dx;
  94       if (fabs(xtemp) >  xmax) xmax = fabs(xtemp);
  95       if (method == 0) xnew[i] = xtemp;
  96       if (method == 1) x[i] = xtemp;
  97     }
  98     if (method == 0)
  99       for (i=0; i<n; i++) x[i] = xnew[i];
 100     (*niter)++;
 101     error /= xmax;
 102     printf("(NITER,ERROR)=%d %g\n", *niter, error);
 103     if (*niter == m)
 104       return;
 105   } while (error > eps);
 106 }