A..0.0.13 chapter6/p132.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 132
   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 double rnd(int *i);
  25 void trid(double a[], double b[], double c[], double f[], int n, int id);
  26 void inpr(double x[], double y[], int n, double *s);
  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   a[n-1] = - 2.0 / h2;
  54   b[n-1] = 2.0 / h2 - shift;
  55   // CALCULATION OF EIGENPAIR
  56   shinv(a,b,c,u,v,eps,h,shift,n,m,&eval,&niter,&npivot);
  57   // OUTPUT
  58   printf("NO. OF ITERATIONS=%3d   EPS=%12.4e   SHIFT=%12.4e   NPIVOT=%3d\n",
  59          niter, eps, shift, npivot);
  60   printf("EIGENVALUE=%12.4e\n", eval);
  61   printf("** NODAL VALUES OF EIGENFUNCTION **\n");
  62   for (i = 0; i < 5; i++)
  63     printf("  I     U(I)    ");
  64   printf("\n");
  65   for (i=0; i<n; i++) {
  66     printf("%3d%12.4e ", i+1, u[i]);
  67     if ((i+1) % 5 == 0 || (i+1) == n)
  68       printf("\n");
  69   }
  70 }
  71 
  72 void shinv(double a[], double b[], double c[],
  73            double x[], double y[],
  74            double eps, double h, double shift,
  75            int n, int m, double *eval, int *niter, int *npivot)
  76 {
  77   int j, i;
  78   double s, t, enew, error;
  79   // STARTING QUANTITIES
  80   j=1;
  81   for (i=0; i<n; i++)
  82     x[i] = rnd(&j);
  83   *niter=0;
  84   *eval=shift;
  85   // ITERATION
  86   do {
  87     for (i=0; i<n; i++)
  88       y[i] = x[i];
  89     trid(a,b,c,y,n,*niter);
  90     inpr(x,y,n,&s);
  91     inpr(y,y,n,&t);
  92     enew=s/t+shift;
  93     t=sqrt(h*t);
  94     for (i=0; i<n; i++)
  95       x[i] = y[i] / t;
  96     error=fabs((enew-*eval)/enew);
  97     *eval=enew;
  98     (*niter)++;
  99     printf("(NITER,EVAL,ERROR)=%d %g %g\n", *niter, *eval, error);
 100     if (*niter == m) break;
 101   }
 102   while (error > eps);
 103   *npivot=0;
 104   for (i=0; i<n; i++)
 105     if (b[i] < 0.0) (*npivot)++;
 106 }
 107 
 108 
 109 double rnd(int *i)
 110 {
 111   // PSEUDO-RANDOM NUMBERS
 112   *i = 12869 * (*i) + 6925;
 113   *i = *i % (1<<15);
 114   return  (double) *i / (1 << 15);
 115 }
 116 
 117 void trid(double a[], double b[], double c[], double f[], int n, int id)
 118 {
 119   int i;
 120   double aa;
 121   // FORWARD ELIMINATION
 122   for (i=0; i<n-1; i++) {
 123     aa = a[i+1] / b[i];
 124     if (id == 0) b[i+1] -= aa * c[i];
 125     f[i+1] -= aa * f[i];
 126   }
 127   // BACKWARD SUBSTITUTION
 128   f[n-1] /= b[n-1];
 129   for (i=n-2; i>=0; i--)
 130     f[i] = (f[i] - c[i]*f[i+1]) / b[i];
 131 }
 132 
 133 void inpr(double x[], double y[], int n, double *s)
 134 {
 135   int i;
 136   // INNER PRODUCT
 137   *s=0.0;
 138   for (i=0; i<n-1; i++)
 139     *s += x[i] * y[i];
 140   *s += x[n-1] * y[n-1] / 2.0;
 141 }