A..0.0.7 chapter3/p73.c

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