C FIGURE 5.11.1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) EXTERNAL FUNA,FUNB CHARACTER*1 HORS C N_ = NUMBER OF SUBINTERVALS C M_ = NUMBER OF UNKNOWNS PARAMETER (N_=4,M_=2*N_+1) C BREAK POINTS ARE XPTS(K), K=0,N COMMON /BLKX/ N,XPTS(0:N_) C UNKNOWNS ARE A(I), I=1,M COMMON /BLKA/ M,A(M_) DIMENSION AMAT(M_,-3:3),F(M_),AOLD(M_) EXACT(X) = EXP(X)*SIN(2.0287579*X) N = N_ M = M_ ALPHA = 4.0 DO 5 K=0,N XPTS(K) = DBLE(K)/N 5 CONTINUE C INITIAL VALUES ASSIGNED DO 10 I=1,M CALL LOCATE(I,HORS,K) IF (HORS.EQ.'H') A(I) = 1.0 IF (HORS.EQ.'S') A(I) = 0.0 10 CONTINUE C DO ITMAX INVERSE POWER ITERATIONS ITMAX = 15 DO 45 ITER=1,ITMAX C SAVE PREVIOUS A VECTOR DO 15 I=1,M AOLD(I) = A(I) 15 CONTINUE C CALCULATE MATRIX A-ALPHA*B, AND F=B*AOLD L = 3 DO 25 K=1,M F(K) = 0.0 DO 20 I=MAX(K-L,1),MIN(K+L,M) AMAT(K,I-K) = GAUSS(FUNA,K,I) BMAT = GAUSS(FUNB,K,I) F(K) = F(K) + BMAT*A(I) AMAT(K,I-K) = AMAT(K,I-K) - ALPHA*BMAT 20 CONTINUE 25 CONTINUE C SOLVE LINEAR SYSTEM USING BAND SOLVER CALL LBAND(AMAT,A,F,M,L) C ESTIMATE EIGENVALUE ANUM = 0.0 ADEN = 0.0 DO 30 I=1,M ANUM = ANUM + F(I)*AOLD(I) ADEN = ADEN + F(I)*A(I) 30 CONTINUE EIGEN = ANUM/ADEN + ALPHA C CALCULATE MAXIMUM ERROR ERMAX = 0.0 DO 35 J=0,100 X = J/100.D0 ERR = & ABS(USOL(0,X)/USOL(0,1.D0) - EXACT(X)/EXACT(1.D0)) ERMAX = MAX(ERMAX,ERR) 35 CONTINUE PRINT 40, ITER,EIGEN,ERMAX 40 FORMAT (I5,2E17.8) 45 CONTINUE STOP END FUNCTION FUNA(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C INTEGRAND FOR (K,I)-TH ELEMENT OF C "A" MATRIX PK = PHI(K,0,X) PKx = PHI(K,1,X) PI = PHI(I,0,X) PIx = PHI(I,1,X) FUNA = EXP(-2.0*X)*(PKx*PIx-PK*PI) RETURN END FUNCTION FUNB(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C INTEGRAND FOR (K,I)-TH ELEMENT OF C "B" MATRIX PK = PHI(K,0,X) PI = PHI(I,0,X) FUNB = EXP(-2.0*X)*PK*PI RETURN END SUBROUTINE LOCATE(I,HORS,K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 HORS C DETERMINE THE TYPE ('H' OR 'S') OF BASIS C FUNCTION I, AND THE POINT (XPTS(K)) ON C WHICH IT IS CENTERED. IF (MOD(I,2).EQ.0) HORS = 'H' IF (MOD(I,2).EQ.1) HORS = 'S' K = I/2 RETURN END FUNCTION OMEGA(IDER,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C EVALUATE OMEGA (IDER=0) OR ITS IDER-TH C DERIVATIVE AT X OMEGA = 0.0 RETURN END