C FIGURE 5.3.3a IMPLICIT DOUBLE PRECISION(A-H,O-Z) EXTERNAL FUNJ,FUNF CHARACTER*1 HORS C N_ = NUMBER OF SUBINTERVALS C M_ = NUMBER OF UNKNOWNS PARAMETER (N_=4,M_=2*N_) 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 AJACOB(M_,-3:3),F(M_),DELTA(M_) DATA PI/3.14159265358979312D0/ N = N_ M = M_ DO 5 K=0,N XPTS(K) = K*PI/N 5 CONTINUE C INITIAL GUESS FOR NEWTON METHOD DO 10 I=1,M CALL LOCATE(I,HORS,K) IF (HORS.EQ.'H') A(I) = 0.5 - XPTS(K)/PI IF (HORS.EQ.'S') A(I) = -1.0/PI 10 CONTINUE C DO ITMAX NEWTON ITERATIONS ITMAX = 15 DO 40 ITER=1,ITMAX C CALCULATE F VECTOR AND JACOBIAN L = 3 DO 20 K=1,M F(K) = GAUSS(FUNF,K,0) DO 15 I=MAX(K-L,1),MIN(K+L,M) AJACOB(K,I-K) = GAUSS(FUNJ,K,I) 15 CONTINUE 20 CONTINUE C SOLVE LINEAR SYSTEM USING BAND SOLVER CALL LBAND(AJACOB,DELTA,F,M,L) C UPDATE SOLUTION VECTOR DO 25 I=1,M A(I) = A(I)-DELTA(I) 25 CONTINUE C CALCULATE MAXIMUM ERROR ERMAX = 0.0 DO 30 J=0,100 X = J*PI/100.D0 ERR = ABS(USOL(0,X) - SIN(X+PI/6.D0)) ERMAX = MAX(ERMAX,ERR) 30 CONTINUE PRINT 35, ITER,ERMAX 35 FORMAT (I5,E15.5) 40 CONTINUE STOP END FUNCTION FUNF(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C INTEGRAND FOR K-TH ELEMENT OF F VECTOR C (I IS NOT USED) U = USOL(0,X) Ux = USOL(1,X) PK = PHI(K,0,X) PKx = PHI(K,1,X) FUNF = -Ux*PKx -Ux**2*PK - U**2*PK + U*PK + PK RETURN END FUNCTION FUNJ(X,K,I) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C INTEGRAND FOR (K,I)-TH ELEMENT C OF JACOBIAN MATRIX U = USOL(0,X) Ux = USOL(1,X) PK = PHI(K,0,X) PKx = PHI(K,1,X) PI = PHI(I,0,X) PIx = PHI(I,1,X) FUNJ = -PKx*PIx - 2*Ux*PK*PIx - 2*U*PK*PI + 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 = 'S' IF (MOD(I,2).EQ.1) HORS = 'H' IF (I.EQ.1) HORS = 'S' K = I/2 RETURN END FUNCTION OMEGA(IDER,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /BLKX/ N,XPTS(0:1) C EVALUATE OMEGA (IDER=0) OR ITS IDER-TH C DERIVATIVE AT X OMEGA = 0.5*HERM('H',0,IDER,X) - 0.5*HERM('H',N,IDER,X) RETURN END