A..0.0.12 chapter6/p125.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 125
   2 /* **        EIGEN ANALYSIS OF       **
   3    ** ORDINARY DIFFERENTIAL EQUATION **
   4    N : NO．OF DIVISIONS
   5    SHIFT : SHIFT PARAMETER
   6    EPS : THRESHOLD VALUE FOR ERROR
   7    M : MAXIMUM NO. OF ITERATIONS
   8    EVAL : APPROXIMATE EIGENVALUE
   9    H : MESH SIZE
  10    A,B,C : COEFFICIENT MATRIX
  11    U : EIGENVECTOR
  12    NITER : NO. OF ITERATIONS
  13    NPIVOT : NO. OF NEGATIVE PIVOTS
  14 */
  15 
  16 #include <stdio.h>
  17 #include <stdlib.h>
  18 #include <math.h>
  19 
  20 void shinv(double a[], double b[], double c[],
  21        double x[], double y[],
  22        double eps, double h, double shift,
  23        int n, int m, double *eval, int *niter, int *npivot);
  24 void trid(double a[], double b[], double c[], double f[], int n, int id);
  25 void inpr(double x[], double y[], int n, double *s);
  26 double rnd(int *i);
  27 
  28 #define NDIM (100)
  29 
  30 int main(void)
  31 {
  32   int n, m, i, niter, npivot;
  33   double shift, eps, h, h2, eval;
  34   double a[NDIM], b[NDIM], c[NDIM], u[NDIM], v[NDIM];
  35   // INPUT
  36   printf("INPUT : N=\n");
  37   scanf("%d", &n);
  38   if (n <= 1) exit(1);
  39   printf("INPUT : SHIFT=\n");
  40   scanf("%lf", &shift);
  41   printf("INPUT : EPS=\n");
  42   scanf("%lf", &eps);
  43   printf("INPUT : M=\n");
  44   scanf("%d", &m);
  45   // CALCULATION OF COEFFICIENT MATRIX
  46   h = 1.0 / n;
  47   h2 = h * h;
  48   for (i=0; i<n-1; i++) {
  49     a[i] = -1.0 / h2;
  50     b[i] = 2.0 / h2 - shift;
  51     c[i] = -1.0 / h2;
  52   }
  53   // CALCULATION OF EIGENPAIR
  54   shinv(a,b,c,u,v,eps,h,shift,n-1,m,&eval,&niter,&npivot);
  55   // OUTPUT
  56   printf("NO. OF ITERATIONS=%3d   EPS=%12.4e   SHIFT=%12.4e   NPIVOT=%3d\n",
  57      niter, eps, shift, npivot);
  58   printf("EIGENVALUE=%12.4e\n", eval);
  59   printf("** NODAL VALUES OF EIGENFUNCTION **\n");
  60   for (i = 0; i < 5; i++)
  61     printf("  I     U(I)    ");
  62   printf("\n");
  63   for (i=0; i<n-1; i++) {
  64     printf("%3d%12.4e ", i+1, u[i]);
  65     if ((i+1) % 5 == 0 || (i+1) == n-1)
  66       printf("\n");
  67   }
  68 }
  69 
  70 void shinv(double a[], double b[], double c[],
  71        double x[], double y[],
  72        double eps, double h, double shift,
  73        int n, int m, double *eval, int *niter, int *npivot)
  74 {
  75   int j, i;
  76   double s, t, enew, error;
  77   // STARTING QUANTITIES
  78   j=1;
  79   for (i=0; i<n; i++)
  80     x[i] = rnd(&j);
  81   *niter=0;
  82   *eval=shift;
  83   // ITERATION
  84   do {
  85     for (i=0; i<n; i++)
  86       y[i] = x[i];
  87     trid(a,b,c,y,n,*niter);
  88     inpr(x,y,n,&s);
  89     inpr(y,y,n,&t);
  90     enew=s/t+shift;
  91     t=sqrt(h*t);
  92     for (i=0; i<n; i++)
  93       x[i] = y[i] / t;
  94     error=fabs((enew-*eval)/enew);
  95     *eval=enew;
  96     (*niter)++;
  97     printf("(NITER,EVAL,ERROR)=%d %g %g\n", *niter, *eval, error);
  98     if (*niter == m) break;
  99   }
 100   while (error > eps);
 101   *npivot=0;
 102   for (i=0; i<n; i++)
 103     if (b[i] < 0.0) (*npivot)++;
 104 }
 105 
 106 double rnd(int *i)
 107 {
 108   // PSEUDO-RANDOM NUMBERS
 109   *i = 12869 * (*i) + 6925;
 110   *i = *i % (1<<15);
 111   return  (double) *i / (1 << 15);
 112 }
 113 
 114 void trid(double a[], double b[], double c[], double f[], int n, int id)
 115 {
 116   int i;
 117   double aa;
 118   // FORWARD ELIMINATION
 119   for (i=0; i<n-1; i++) {
 120     aa = a[i+1] / b[i];
 121     if (id == 0) b[i+1] -= aa * c[i];
 122     f[i+1] -= aa * f[i];
 123   }
 124   // BACKWARD SUBSTITUTION
 125   f[n-1] /= b[n-1];
 126   for (i=n-2; i>=0; i--)
 127     f[i] = (f[i] - c[i]*f[i+1]) / b[i];
 128 }
 129 
 130 void inpr(double x[], double y[], int n, double *s)
 131 {
 132   int i;
 133   // INNER PRODUCT
 134   *s=0.0;
 135   for (i=0; i<n; i++)
 136     *s += x[i] * y[i];
 137 }