C FIGURE 5.8.1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) EXTERNAL FUNF,FUNA,FUNB CHARACTER*1 HORS C N_ = NUMBER OF SUBINTERVALS C M_ = NUMBER OF UNKNOWNS C NSTEP = NUMBER OF TIME STEPS PARAMETER (N_=10,M_=2*N_+1,NSTEP=200) 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_) C PASS T(k+1) THROUGH COMMON COMMON /BLKT/ TKP1 DIMENSION AMAT(M_,-3:3),F(M_) N = N_ M = M_ DO 5 K=0,N XPTS(K) = DBLE(K)/N 5 CONTINUE C INITIAL VALUES SET DO 10 I=1,M CALL LOCATE(I,HORS,K) IF (HORS.EQ.'H') A(I) = EXP(-XPTS(K)) IF (HORS.EQ.'S') A(I) = -EXP(-XPTS(K)) 10 CONTINUE DT = 1.D0/NSTEP DO 40 KT=1,NSTEP TKP1 = KT*DT C BACKWARD DIFFERENCE METHOD REQUIRES C SOLUTION OF (BMAT + DT*AMAT)*a(k+1) C = DT*b + BMAT*a(k) EACH TIME STEP L = 3 DO 20 K=1,M F(K) = DT*GAUSS(FUNF,K,0) DO 15 I=MAX(K-L,1),MIN(K+L,M) AMATKI = GAUSS(FUNA,K,I) + & EXP(0.5*TKP1)*PHI(K,0,1.D0)*PHI(I,0,1.D0) BMATKI = GAUSS(FUNB,K,I) F(K) = F(K) + BMATKI*A(I) AMAT(K,I-K) = BMATKI + DT*AMATKI 15 CONTINUE 20 CONTINUE C SOLVE LINEAR SYSTEM USING BAND SOLVER CALL LBAND(AMAT,A,F,M,L) C CALCULATE MAXIMUM ERROR ERMAX = 0.0 DO 30 J=0,100 X = J/100.D0 ERR = ABS(USOL(0,X) - EXP(TKP1-X)) ERMAX = MAX(ERMAX,ERR) 30 CONTINUE PRINT 35, TKP1,ERMAX 35 FORMAT (F8.3,E15.5) 40 CONTINUE STOP END FUNCTION FUNF(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /BLKT/ T C INTEGRAND FOR K-TH ELEMENT OF "b" VECTOR C (I IS NOT USED) OM = OMEGA(0,X) OMt = OMEGA(0,X) OMx = OMEGA(1,X) PK = PHI(K,0,X) PKx = PHI(K,1,X) FUNF = EXP(0.5*X**2*T)*(-OMt*PK - OMx*PKx + X*T*OM*PK) RETURN END FUNCTION FUNA(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /BLKT/ T C INTEGRAND FOR (K,I)-TH ELEMENT C OF "A" MATRIX PK = PHI(K,0,X) PKx = PHI(K,1,X) PI = PHI(I,0,X) PIx = PHI(I,1,X) FUNA = EXP(0.5*X**2*T)*(PKx*PIx - X*T*PK*PI) RETURN END FUNCTION FUNB(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /BLKT/ T C INTEGRAND FOR (K,I)-TH ELEMENT C OF "B" MATRIX PK = PHI(K,0,X) PI = PHI(I,0,X) FUNB = EXP(0.5*X**2*T)*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) COMMON /BLKT/ T C EVALUATE OMEGA (IDER=0) OR ITS IDER-TH C DERIVATIVE AT X OMEGA = EXP(T)*HERM('H',0,IDER,X) RETURN END