C FIGURE 5.3.3b FUNCTION USOL(IDER,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /BLKA/ M,A(1) C EVALUATE THE APPROXIMATE SOLUTION (IDER=0) C OR ITS IDER-TH DERIVATIVE AT X USOL = OMEGA(IDER,X) DO 10 I=1,M USOL = USOL + A(I)*PHI(I,IDER,X) 10 CONTINUE RETURN END FUNCTION PHI(I,IDER,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 HORS C EVALUATE THE I-TH BASIS FUNCTION (IDER=0) C OR ITS IDER-TH DERIVATIVE AT X CALL LOCATE(I,HORS,K) PHI = HERM(HORS,K,IDER,X) RETURN END FUNCTION HERM(HORS,K,IDER,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 HORS COMMON /BLKX/ N,XPTS(0:1) C EVALUATE AT X THE IDER-TH DERIVATIVE OF C THE CUBIC HERMITE "H" OR "S" FUNCTION C CENTERED AT XPTS(K). XK = XPTS(K) HH = XPTS(N)-XPTS(0) IF (K.GT.0) XKM1 = XPTS(K-1) IF (K.LE.0) XKM1 = XPTS(0)-HH IF (K.LT.N) XKP1 = XPTS(K+1) IF (K.GE.N) XKP1 = XPTS(N)+HH HERM = 0.0 IF (X.LT.XKM1) RETURN IF (X.GT.XKP1) RETURN IF (X.LE.XK) THEN H = XK-XKM1 S = (X-XKM1)/H IF (HORS.EQ.'H') THEN IF (IDER.EQ.0) HERM = 3.0*S**2 - 2.0*S**3 IF (IDER.EQ.1) HERM = 6.0*S/H - 6.0*S**2/H IF (IDER.EQ.2) HERM = 6.0/H**2 - 12.0*S/H**2 ELSE IF (IDER.EQ.0) HERM = -H*S**2 + H*S**3 IF (IDER.EQ.1) HERM = -2.0*S + 3.0*S**2 IF (IDER.EQ.2) HERM = -2.0/H + 6.0*S/H ENDIF ELSE H = XKP1-XK S = (XKP1-X)/H IF (HORS.EQ.'H') THEN IF (IDER.EQ.0) HERM = 3.0*S**2 - 2.0*S**3 IF (IDER.EQ.1) HERM = -6.0*S/H + 6.0*S**2/H IF (IDER.EQ.2) HERM = 6.0/H**2 - 12.0*S/H**2 ELSE IF (IDER.EQ.0) HERM = H*S**2 - H*S**3 IF (IDER.EQ.1) HERM = -2.0*S + 3.0*S**2 IF (IDER.EQ.2) HERM = 2.0/H - 6.0*S/H ENDIF ENDIF RETURN END FUNCTION GAUSS(F,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 HORS DIMENSION S(3),W(3) COMMON /BLKX/ N,XPTS(0:1) C USE THE GAUSS 3 POINT FORMULA TO C INTEGRATE F(X,K,I). DATA S & /0.1127016653792583D0,0.5000000000000000D0,0.8872983346207417D0/ DATA W & /0.2777777777777778D0,0.4444444444444444D0,0.2777777777777778D0/ GAUSS = 0.0 C THE FUNCTION F INVOLVES PHI(K,...), WHICH C IS NONZERO ONLY FROM XPTS(IK-1) TO C XPTS(IK+1) CALL LOCATE(K,HORS,IK) DO 10 J=MAX(IK,1),MIN(IK+1,N) H = XPTS(J)-XPTS(J-1) DO 5 JX=1,3 X = XPTS(J-1)+S(JX)*H GAUSS = GAUSS + H*W(JX)*F(X,K,I) 5 CONTINUE 10 CONTINUE RETURN END