PROGRAM FIXPOINT C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.1 (Fixed Point Iteration). C Section 2.1, Iteration for Solving x = g(x), Page 51 C PARAMETER(Tol=0.00000001,Max=500) INTEGER Cond,K REAL P0,Pnew,Pterm EXTERNAL G 10 CALL GETPOINT(P0,Pterm) CALL ITERATE(G,Pterm,Max,Tol,Pnew,Cond,K) CALL RESULTS(P0,Max,Tol,Pnew,Pterm,Cond,K) WRITE(9,*) ' WANT TO DO ANOTHER ONE ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END REAL FUNCTION G(X) REAL X G=0.9 + X - 0.4*X*X RETURN END SUBROUTINE PRINTFN WRITE(9,*)'G(X) = 0.9 + X - 0.4*X*X' RETURN END SUBROUTINE ITERATE(G,Pterm,Max,Tol,Pnew,Cond,K) PARAMETER(Big=1E10,Small=1E-20) INTEGER Cond,K,Max REAL Pnew,Pterm,Tol REAL Dx,Dg,Pold,RelErr,Slope K=0 RelErr=1 Pnew=G(Pterm) WHILE (RelErr.GE.Tol).AND.(K.LE.Max) &.AND.(ABS(Pnew).LT.Big) Pold=Pterm Pterm=Pnew Pnew=G(Pterm) Dg=Pnew-Pterm RelErr=ABS(Dg)/(ABS(Pnew)+Small) K=K+1 WRITE(9,1000) K,Pnew REPEAT IF (Dg.EQ.0) THEN Slope=0 ELSE Dx=Pterm-Pold IF (Dx.NE.0) THEN Slope=Dg/Dx ELSE Slope=6.023E23 ENDIF ENDIF IF (ABS(Slope).LT.1) THEN Cond=1 IF (Slope.LT.0) Cond=2 ELSE Cond=3 IF (Slope.LT.0) Cond=4 ENDIF IF (RelErr.LT.Tol) THEN IF ((Cond.EQ.3).OR.(Cond.EQ.4)) Cond=5 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7) END SUBROUTINE XITERATE(G,Pterm,Max,Tol,Pnew,Cond,K) C This subroutine uses simulated WHILE loop(s). PARAMETER(Big=1E10,Small=1E-20) INTEGER Cond,K,Max REAL Pnew,Pterm,Tol REAL Dx,Dg,Pold,RelErr,Slope K=0 RelErr=1 Pnew=G(Pterm) 10 IF ((RelErr.GE.Tol).AND.(K.LE.Max)) THEN Pold=Pterm Pterm=Pnew Pnew=G(Pterm) Dg=Pnew-Pterm RelErr=ABS(Dg)/(ABS(Pnew)+Small) K=K+1 IF ((Pnew.LT.(-Big)).OR.(Big.LT.Pnew)) GOTO 20 WRITE(9,1000) K,Pnew GOTO 10 ENDIF 20 CONTINUE !Precludes Overflow IF (Dg.EQ.0) THEN Slope=0 ELSE Dx=Pterm-Pold IF (Dx.NE.0) THEN Slope=Dg/Dx ELSE Slope=6.023E23 ENDIF ENDIF IF (ABS(Slope).LT.1) THEN Cond=1 IF (Slope.LT.0) Cond=2 ELSE Cond=3 IF (Slope.LT.0) Cond=4 ENDIF IF (RelErr.LT.Tol) THEN IF ((Cond.EQ.3).OR.(Cond.EQ.4)) Cond=5 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7) END SUBROUTINE GETPOINT(P0,Pterm) REAL P0,Pterm DO 10 I=1,17 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'FIXED POINT ITERATION IS USED TO FIND' WRITE(9,*)' ' WRITE(9,*)'A FIXED POINT OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'ENTER THE STARTING VALUE P0 = ' READ(9,*) P0 Pterm=P0 WRITE(9,*)' ' RETURN END SUBROUTINE RESULTS(P0,Max,Tol,Pnew,Pterm,Cond,K) PARAMETER(Small=1E-20) INTEGER Cond,I,K,Max REAL Delta,P0,Pnew,Pterm,RelErr,Tol DO 10 I=1,12 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'FIXED POINT ITERATION WAS USED TO FIND A FIXED POINT OF +' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'THE STARTING POINT WAS P(0) =',P0 WRITE(9,*)' ' Delta=ABS(Pnew-Pterm) RelErr=Delta/(ABS(Pnew)+Small) IF (RelErr.LE.Tol) THEN WRITE(9,*)'AFTER ',K,' ITERATIONS THE FIXED POINT WAS FOUND' WRITE(9,*)' ' WRITE(9,*)' P =',Pnew WRITE(9,*)' ' WRITE(9,*)'ITERATION WAS SUCCESSFUL.' ELSE WRITE(9,*)'THE LOCATION OF THE FIXED POINT IS IN DOUBT.' WRITE(9,*)' ' WRITE(9,*)'THE APPROXIMATION AFTER ',K,' ITERATIONS IS' WRITE(9,*)' ' WRITE(9,*)'P(',K,') =',Pnew ENDIF WRITE(9,*)' ' WRITE(9,*)' ',Pnew,' = G(',Pterm,' )' WRITE(9,*)' ' WRITE(9,*)'CONSECUTIVE ITERATES ARE WITHIN',Delta WRITE(9,*)' ' IF (K.GE.Max) THEN WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED' WRITE(9,*)' ' ENDIF IF (Cond.EQ.1) THEN WRITE(9,*)'THE SEQUENCE EXHIBITS MONOTONE CONVERGENCE.' ELSEIF (Cond.EQ.2) THEN WRITE(9,*)'THE SEQUENCE EXHIBITS OSCILLATING CONVERGENCE.' ELSEIF (Cond.EQ.3) THEN WRITE(9,*)'THE SEQUENCE EXHIBITS MONOTONE DIVERGENCE.' ELSEIF (Cond.EQ.4) THEN WRITE(9,*)'THE SEQUENCE EXHIBITS OSCILLATING DIVERGENCE.' ENDIF DO 20 I=1,1 WRITE(9,*)' ' 20 CONTINUE RETURN END PROGRAM BISECTION C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.2 (Bisection Method). C Section 2.2, Bracketing Methods for Locating a Root, Page 61 C PARAMETER(Delta=5E-7) INTEGER Cond,KL,KR REAL A,A1,B,B1,C,D CHARACTER ANS*1 EXTERNAL F 10 CALL GETPTS(A,B) A1=A B1=B CALL BISEC(F,A1,B1,Delta,C,D,Cond,KL,KR) CALL RESULT(F,A,B,Delta,C,D,Cond,KL,KR) WRITE(9,*)' ' WRITE(9,*)'WANT TO TRY ANOTHER INTERVAL? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END REAL FUNCTION F(X) F = SIN(X)/COS(X) RETURN END SUBROUTINE PRINTFN WRITE(9,*)'F(X) = SIN(X)/COS(X)' RETURN END SUBROUTINE BISEC(F,A,B,Delta,C,D,Cond,KL,KR) PARAMETER(Big=1E5) INTEGER Cond,K,KL,KR,Max REAL A,B,C,D,Delta,YA,YB,YC EXTERNAL F K=0 KL=0 KR=0 YA=F(A) YB=F(B) D=B-A Cond=0 Max=1+INT((ALOG(D)-ALOG(Delta))/ALOG(2.0)) IF (YA*YB.GE.0) THEN Cond=1 RETURN ENDIF WHILE (Cond.EQ.0).AND.(K.LT.Max) C=(A+B)/2.0 YC=F(C) IF (YC.EQ.0) THEN A=C B=C Cond=2 ELSEIF ((YB*YC).GT.0) THEN B=C YB=YC KR=KR+1 ELSE A=C YA=YC KL=KL+1 ENDIF K=K+1 WRITE(9,1000) K,A,C,B REPEAT D=B-A IF (D.LT.Delta) THEN IF (Cond.NE.2) Cond=3 IF ((ABS(YA).GT.Big).AND.(ABS(YB).GT.Big)) Cond=4 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7,4X,F15.7) END SUBROUTINE XBISEC(F,A,B,Delta,C,D,Cond,KL,KR) C This subroutine uses simulated WHILE loop(s). PARAMETER(Big=1E5) INTEGER Cond,K,KL,KR,Max REAL A,B,C,D,Delta,YA,YB,YC EXTERNAL F K=0 KL=0 KR=0 YA=F(A) YB=F(B) D=B-A Cond=0 Max=1+INT((ALOG(D)-ALOG(Delta))/ALOG(2.0)) IF (YA*YB.GE.0) THEN Cond=1 RETURN ENDIF 10 IF ((Cond.EQ.0).AND.(K.LT.Max)) THEN C=(A+B)/2.0 YC=F(C) IF (YC.EQ.0) THEN A=C B=C Cond=2 ELSEIF ((YB*YC).GT.0) THEN B=C YB=YC KR=KR+1 ELSE A=C YA=YC KL=KL+1 ENDIF K=K+1 WRITE(9,1000) K,A,C,B GOTO 10 ENDIF D=B-A IF (D.LT.Delta) THEN IF (Cond.NE.2) Cond=3 IF ((ABS(YA).GT.Big).AND.(ABS(YB).GT.Big)) Cond=4 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7,4X,F15.7) END SUBROUTINE GETPTS(A,B) INTEGER I REAL A,B DO 10 I=1,17 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'BISECTION METHOD IS USED TO FIND A ZERO THE FUNCTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'THAT LIES IN THE INTERVAL [A,B]' WRITE(9,*)' ' WRITE(9,*)'ENTER LEFT ENDPOINT A = ' READ(9,*) A WRITE(9,*)'ENTER RIGHT ENDPOINT B = ' READ(9,*) B WRITE(9,*)' ' RETURN END SUBROUTINE RESULT(F,A,B,Delta,C,D,Cond,KL,KR) INTEGER Cond,I,J,K,KL,KR REAL A,B,Delta,C,D EXTERNAL F DO 10 I=1,16 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE BISECTION METHOD WAS USED TO FIND A ZERO OF THE FUN +CTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'IN THE INTERVAL [',A,' ,',B,' ]' WRITE(9,*)' ' K=KL+KR IF (Cond.EQ.1) THEN WRITE(9,*)'THE VALUES F(A) AND F(B)' WRITE(9,*)'DO NOT DIFFER IN SIGN.' WRITE(9,*)' ' WRITE(9,*)'F(',A,' ) =',F(A) WRITE(9,*)'F(',B,' ) =',F(B) WRITE(9,*)' ' WRITE(9,*)'THE BISECTION METHOD CANNOT BE USED.' ELSEIF ((Cond.EQ.2).OR.(Cond.EQ.3)) THEN WRITE(9,*)'THE BISECTION METHOD TOOK ',K,' ITERATIONS.' WRITE(9,*)' ' WRITE(9,*)'THERE WERE',KR,' SQUEEZES FROM THE RIGHT' WRITE(9,*)' ' WRITE(9,*)' AND',KL,' SQUEEZES FROM THE LEFT.' WRITE(9,*)' ' WRITE(9,*)' THE ROOT IS C =',C WRITE(9,*)' ' WRITE(9,*)'THE ACCURACY IS +-',D WRITE(9,*)' ' WRITE(9,*)'FUNCTION VALUE F(',C,' ) =',F(C) WRITE(9,*)' ' IF (Cond.EQ.2) THEN WRITE(9,*)'THE FUNCTION VALUE IS EXACTLY ZERO!' ENDIF ELSEIF ((Cond.EQ.0).OR.(Cond.EQ.4)) THEN WRITE(9,*)'THE CURRENT ITERATE IS C(',K,') =',C WRITE(9,*)' ' WRITE(9,*)'THE CURRENT INTERVAL WIDTH IS ',D WRITE(9,*)' ' WRITE(9,*)'THE CURRENT FUNCTION VALUE F(',C,' ) =',F(C) WRITE(9,*)' ' IF (Cond.EQ.0) THEN WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ENDIF IF (Cond.EQ.4) THEN WRITE(9,*)'SURPRISE, A POLE WAS FOUND INSTEAD OF A ZERO!' ENDIF ENDIF WRITE(9,*)' ' RETURN END PROGRAM FALSEPOS C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.3 (False position or Regula Falsi Method). C Section 2.2, Bracketing Methods for Locating a Root, Page 62 C PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=99) INTEGER Cond,K REAL A,A1,B,B1,C,DX,YC CHARACTER ANS*1 EXTERNAL F 10 CALL GETPTS(A,B) A1=A B1=B CALL REGULA(F,A1,B1,Delta,Epsilon,Max,C,YC,D,K,Cond) CALL RESULT(F,A,B,Delta,Epsilon,Max,C,YC,D,K,Cond) WRITE(9,*)' ' WRITE(9,*)'WANT TO TRY ANOTHER INTERVAL? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END REAL FUNCTION F(X) F = SIN(X)/COS(X) RETURN END SUBROUTINE PRINTFN WRITE(9,*)'F(X) = SIN(X)/COS(X)' RETURN END SUBROUTINE REGULA(F,A,B,Delta,Epsilon,Max,C,YC,D,K,Cond) PARAMETER(Big=1E3) INTEGER Cond,K,Max REAL A,B,C,D,Delta,DX,Epsilon,YC EXTERNAL F K=0 YA=F(A) YB=F(B) D=B-A Cond=0 IF (YA*YB.GE.0) THEN Cond=1 RETURN ENDIF WHILE ((Cond.EQ.0).AND.(K.LT.Max)) DXA=YA*(B-A)/(YB-YA) DXB=YB*(B-A)/(YB-YA) IF ABS(DXA).LT.ABS(DXB) THEN C=A-DXA D=ABS(DXA) ELSE C=B-DXB D=ABS(DXB) ENDIF YC=F(C) IF (D.LT.Delta) Cond=3 IF (ABS(YC).LT.Epsilon) THEN A=C B=C Cond=4 ELSEIF ((YB*YC).GT.0) THEN B=C YB=YC ELSE A=C YA=YC ENDIF K=K+1 WRITE(9,1000) K,A,C,B REPEAT IF (YC.EQ.0) Cond=2 IF (D.LT.Delta).AND.(Cond.NE.2) Cond=3 IF (Cond.EQ.3).AND.(ABS(YC).LT.Epsilon) Cond=5 IF (ABS(YA).GT.Big).AND.(ABS(YB).GT.Big) Cond=6 IF (Cond.EQ.0).AND.(ABS(YC).LT.Epsilon) Cond=4 PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7,4X,F15.7) END SUBROUTINE XREGULA(F,A,B,Delta,Epsilon,Max,C,YC,D,K,Cond) C This subroutine uses simulated WHILE loop(s). PARAMETER(Big=1E3) INTEGER Cond,K,Max REAL A,B,C,D,Delta,DX,Epsilon,YC EXTERNAL F K=0 YA=F(A) YB=F(B) D=B-A Cond=0 IF (YA*YB.GE.0) THEN Cond=1 RETURN ENDIF 10 IF ((Cond.EQ.0).AND.(K.LT.Max)) THEN DXA=YA*(B-A)/(YB-YA) DXB=YB*(B-A)/(YB-YA) IF ABS(DXA).LT.ABS(DXB) THEN C=A-DXA D=ABS(DXA) ELSE C=B-DXB D=ABS(DXB) ENDIF YC=F(C) IF (D.LT.Delta) Cond=3 IF (ABS(YC).LT.Epsilon) THEN A=C B=C Cond=4 ELSEIF ((YB*YC).GT.0) THEN B=C YB=YC ELSE A=C YA=YC ENDIF K=K+1 WRITE(9,1000) K,A,C,B,D GOTO 10 ENDIF IF (YC.EQ.0) Cond=2 IF (D.LT.Delta) THEN IF (Cond.NE.2) Cond=3 IF ((Cond.EQ.3).AND.(ABS(YC).LT.Epsilon)) Cond=5 IF ((ABS(YA).GT.Big).AND.(ABS(YB).GT.Big)) Cond=6 ENDIF IF ((Cond.EQ.0).AND.(ABS(YC).LT.Epsilon)) Cond=4 PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7,4X,F15.7,4X,F15.7) END SUBROUTINE GETPTS(A,B) INTEGER I REAL A,B DO 10 I=1,17 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE FALSE POSITION METHOD WILL BE USED' WRITE(9,*)'TO FIND A ZERO THE FUNCTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'THAT LIES IN THE INTERVAL [A,B]' WRITE(9,*)' ' WRITE(9,*)'ENTER LEFT ENDPOINT A = ' READ(9,*) A WRITE(9,*)'ENTER RIGHT ENDPOINT B = ' READ(9,*) B WRITE(9,*)' ' RETURN END SUBROUTINE RESULT(F,A,B,Delta,Epsilon,Max,C,YC,D,K,Cond) INTEGER Cond,I,J,K,Max REAL A,B,C,Delta,D,Epsilon,YC EXTERNAL F DO 10 I=1,16 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE FALSE POSITION METHOD WAS USED TO FIND A ZERO' WRITE(9,*)'OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'IN THE INTERVAL [',A,' ,',B,' ]' WRITE(9,*)' ' IF (Cond.EQ.1) THEN WRITE(9,*)'THE VALUES F(A) AND F(B)' WRITE(9,*)'DO NOT DIFFER IN SIGN.' WRITE(9,*)' ' WRITE(9,*)'F(',A,' ) =',F(A) WRITE(9,*)'F(',B,' ) =',F(B) WRITE(9,*)' ' WRITE(9,*)'THE FALSE POSITION METHOD CANNOT BE USED.' ELSEIF ((Cond.GE.2).AND.(Cond.LE.5)) THEN WRITE(9,*)'THE FALSE POSITION METHOD TOOK ',K,' ITERATIONS.' WRITE(9,*)' ' WRITE(9,*)' THE ROOT IS C =',C WRITE(9,*)' ' WRITE(9,*)'THE ACCURACY IS +-',D WRITE(9,*)' ' WRITE(9,*)'FUNCTION VALUE F(',C,' ) =',F(C) WRITE(9,*)' ' IF (Cond.EQ.2) THEN WRITE(9,*)'THE FUNCTION VALUE IS EXACTLY ZERO!' ENDIF IF ((Cond.EQ.3).OR.(Cond.EQ.5)) THEN WRITE(9,*)'THE ABSCISSA IS WITHIN TOLERANCE' ENDIF IF ((Cond.EQ.4).OR.(Cond.EQ.5)) THEN WRITE(9,*)'THE FUNCTION VALUE IS IS WITHIN TOLERANCE' ENDIF ELSEIF ((Cond.EQ.0).OR.(Cond.EQ.6)) THEN WRITE(9,*)'THE CURRENT ITERATE IS C(',K,') =',C WRITE(9,*)' ' WRITE(9,*)'CONSECUTIVE ITERATES DIFFER BY',D WRITE(9,*)' ' WRITE(9,*)'THE CURRENT FUNCTION VALUE F(',C,' ) =',F(C) WRITE(9,*)' ' IF (Cond.EQ.0) THEN WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ENDIF IF (Cond.EQ.6) THEN WRITE(9,*)'SURPRISE, A POLE WAS FOUND INSTEAD OF A ZERO!' ENDIF ENDIF WRITE(9,*)' ' RETURN END PROGRAM APPROXLO C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.4 (Approximate Location of Roots). C Section 2.3, Initial Approximations and Convergence Criteria, Page 70 C PARAMETER(Epsilon=0.05,MN=201,MR=50) INTEGER N,M REAL A,B,R,X,Y CHARACTER ANS*1 DIMENSION X(0:MN),R(0:MR),Y(0:MN) EXTERNAL F 10 CALL INPUT(A,B,N) CALL LOCATE(F,A,B,N,Epsilon,R,M) CALL RESULTS(F,A,B,N,R,M) WRITE(9,*)' ' WRITE(9,*) 'WANT TO TRY ANOTHER INTERVAL? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y'.OR.ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) REAL X F=SIN(X)+1 RETURN END SUBROUTINE PRINTFN WRITE(9,*)'F(X) = SIN(X)+1' RETURN END SUBROUTINE LOCATE(F,A,B,N,Epsilon,R,M) PARAMETER(MN=201,MR=50) INTEGER K,M,N REAL A,B,Epsilon,H,MaxY,MinY,R,Slope,Small,X,Y DIMENSION X(0:MN),R(0:MR),Y(0:MN) EXTERNAL F H=(B-A)/N DO K=0,N X(K)=A+H*K Y(K)=F(X(K)) ENDDO MaxY=Y(0) DO K=1,N IF (MaxY.LT.Y(K)) MaxY=Y(K) ENDDO MinY=Y(0) DO K=1,N IF (MinY.GT.Y(K)) MinY=Y(K) ENDDO Small=(MaxY-MinY)*Epsilon M=0 Y(N+1)=Y(N) DO K=1,N IF ((Y(K-1)*Y(K)).LE.0) THEN M=M+1 R(M)=(X(K-1)+X(K))/2 ENDIF Slope=(Y(K)-Y(K-1))*(Y(K+1)-Y(K)) IF ((ABS(Y(K)).LT.Small).AND.(Slope.LT.0)) THEN M=M+1 R(M)=X(K) ENDIF ENDDO RETURN END SUBROUTINE XLOCATE(F,A,B,N,Epsilon,R,M) C This subroutine uses simulated WHILE loop(s). PARAMETER(MN=201,MR=50) INTEGER K,M,N REAL A,B,Epsilon,H,MaxY,MinY,R,Slope,Small,X,Y DIMENSION X(0:MN),R(0:MR),Y(0:MN) EXTERNAL F H=(B-A)/N DO 10 K=0,N X(K)=A+H*K Y(K)=F(X(K)) 10 CONTINUE MaxY=Y(0) DO 20 K=1,N IF (MaxY.LT.Y(K)) THEN MaxY=Y(K) ENDIF 20 CONTINUE MinY=Y(0) DO 30 K=1,N IF (MinY.GT.Y(K)) THEN MinY=Y(K) ENDIF 30 CONTINUE Small=(MaxY-MinY)*Epsilon M=0 Y(N+1)=Y(N) DO 40 K=1,N IF ((Y(K-1)*Y(K)).LE.0) THEN M=M+1 R(M)=(X(K-1)+X(K))/2 ENDIF Slope=(Y(K)-Y(K-1))*(Y(K+1)-Y(K)) IF ((ABS(Y(K)).LT.Small).AND.(Slope.LT.0)) THEN M=M+1 R(M)=X(K) ENDIF 40 CONTINUE RETURN END SUBROUTINE INPUT(A,B,N) INTEGER I,N REAL A,B DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'APPROXIMATE LOCATIONS FOR THE ZEROS OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'THE INTERVAL [A,B] IS DIVIDED INTO N SUBINTERVALS' WRITE(9,*)' ' WRITE(9,*)'AND A SEARCH IS PERFORMED TO DETERMINE IF F(X)' WRITE(9,*)' ' WRITE(9,*)'CHANGES SIGN OR IS CLOSE TO ZERO ON EACH SUBINTERVAL.' WRITE(9,*)' ' WRITE(9,*)'ENTER LEFT ENDPOINT A = ' READ(9,*) A WRITE(9,*)' ' WRITE(9,*)'ENTER RIGHT ENDPOINT B = ' READ(9,*) B WRITE(9,*)' ' WRITE(9,*)'NUMBER OF SUBINTERVALS N = ' READ(9,*) N RETURN END SUBROUTINE RESULTS(F,A,B,N,R,M) PARAMETER(MR=50) INTEGER I,K,M,N REAL A,B,R DIMENSION R(0:MR) EXTERNAL F DO 10 I=1,17 WRITE(9,*)' ' 10 CONTINUE H=(B-A)/N WRITE(9,*)'WE SEARCHED TO FIND APPROXIMATE LOCATIONS FOR THE ZEROS + OF' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'IN THE INTERVAL [',A,' ,',B,' ].' WRITE(9,*)' ' WRITE(9,*)'THE SEARCH WAS PERFORMED USING ',N,' SUBINTERVALS' WRITE(9,*)' ' WRITE(9,*)'OF EQUAL WIDTH H =',H WRITE(9,*)' ' WRITE(9,*)' ' IF (M.GE.0) THEN WRITE(9,*)'WE FOUND ',M,' APPROXIMATE LOCATION(S).' WRITE(9,*)' ' WRITE(9,*)'THEIR ABSCISSAS X(K) AND ORDINATES Y(K) ARE:' WRITE(9,*)' ' DO 20 K=1,M WRITE(9,*)'X(',K,') =',R(K),' Y(',K,') =',F(R(K)) 20 CONTINUE ELSE WRITE(9,*)'NO APPROXIMATE LOCATIONS FOR ZEROS OF F(X) WERE FOUND +.' WRITE(9,*)' ' WRITE(9,*)'YOU COULD TRY ANOTHER SEARCH USING A SMALLER STEP H,' WRITE(9,*)' ' WRITE(9,*)'OR YOU COULD TRY A SEARCH IS DIFFERENT INTERVAL [A,B] +.' ENDIF RETURN END PROGRAM NEWTON C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.5 (Newton-Raphson Iteration). C Section 2.4, Newton-Raphson and Secant Methods, Page 84 C PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=100) INTEGER Cond,K REAL Dp,P,P0,P1,Y1,RelErr CHARACTER ANS*1 EXTERNAL F,F1 10 CALL INPUT(P0) P=P0 CALL NEWRAP(F,F1,P,Delta,Epsilon,Max,P1,Dp,Y1,Cond,K) CALL RESULT(P0,P1,Dp,Y1,Cond,K) WRITE(9,*)' ' WRITE(9,*) 'WANT TO USE A DIFFERENT STARTING VALUE? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y'.OR.ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) F=X*X*X-3*X+2 RETURN END REAL FUNCTION F1(X) F1=3*X*X-3 RETURN END SUBROUTINE PRINTFN WRITE(9,*)'F(X) = X*X*X-3*X+2' RETURN END SUBROUTINE NEWRAP(F,F1,P0,Delta,Epsilon,Max,P1,Dp,Y1,Cond,K) PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL Delta,Epsilon,Df,Dp,P0,P1,Y0,Y1,RelErr EXTERNAL F,F1 K=0 Cond=0 Y0=F(P0) P1=P0+1 WHILE ((K.LE.Max).AND.(Cond.EQ.0)) Df=F1(P0) IF (Df.EQ.0) THEN Cond=1 Dp=P1-P0 P1=P0 ELSE Dp=Y0/Df P1=P0-Dp ENDIF Y1=F(P1) RelErr=ABS(Dp)/(ABS(P1)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y1).LT.Epsilon) Cond=3 IF (Cond.EQ.2).AND.(ABS(Y1).LT.Epsilon) Cond=4 P0=P1 Y0=Y1 K=K+1 WRITE(9,1000) K,P1,Y1 REPEAT PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE XNEWRAP(F,F1,P0,Delta,Epsilon,Max,P1,Dp,Y1,Cond,K) C This subroutine uses simulated WHILE loop(s). PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL Delta,Epsilon,Df,Dp,P0,P1,Y0,Y1,RelErr EXTERNAL F,F1 K=0 Cond=0 Y0=F(P0) P1=P0+1 10 IF ((K.LE.Max).AND.(Cond.EQ.0)) THEN Df=F1(P0) IF (Df.EQ.0) THEN Cond=1 Dp=P1-P0 P1=P0 ELSE Dp=Y0/Df P1=P0-Dp ENDIF Y1=F(P1) RelErr=ABS(Dp)/(ABS(P1)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y1).LT.Epsilon) Cond=3 IF ((RelErr.LE.Delta).AND.(ABS(Y1).LT.Epsilon)) Cond=4 P0=P1 Y0=Y1 K=K+1 WRITE(9,1000) K,P1,Y1 GOTO 10 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE INPUT(P0) INTEGER I REAL P0 DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE NEWTON-RAPHSON METHOD IS USED' WRITE(9,*)' ' WRITE(9,*)'TO FIND A ZERO OF THE FUNCTION: ' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'ONE INITIAL APPROXIMATION P0 IS NEEDED.' WRITE(9,*)' ' WRITE(9,*)'ENTER P0 = ' READ(9,*) P0 WRITE(9,*)' ' RETURN END SUBROUTINE RESULT(P0,P1,Dp,Y1,Cond,K) INTEGER Cond,I,K REAL P0,P1,Dp,Y1 DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE NEWTON-RAPHSON METHOD WAS USED TO FIND A ZERO OF' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'STARTING WITH THE APPROXIMATION P0 =',P0 WRITE(9,*)' ' WRITE(9,*)'AFTER ',K,' ITERATIONS AN APPROXIMATION FOR THE ZERO IS +:' WRITE(9,*)' ' WRITE(9,*)' P =',P1 WRITE(9,*)' ' WRITE(9,*)' DP =',ABS(Dp),' IS THE ESTIMATED ACCURACY FOR P.' WRITE(9,*)' ' WRITE(9,*)' F(',P1,' ) =',Y1 WRITE(9,*)' ' IF (Y1.EQ.0) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE IS EXACTLY ZERO!' WRITE(9,*)' ' ENDIF IF (Cond.EQ.0) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE' WRITE(9,*)' ' WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ELSEIF (Cond.EQ.1) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE ' WRITE(9,*)' ' WRITE(9,*)'DIVISION BY ZERO WAS ENCOUNTERED.' ELSEIF (Cond.EQ.2) THEN WRITE(9,*)'THE APPROXIMATION P IS' WRITE(9,*)' ' WRITE(9,*)'WITHIN THE DESIRED TOLERANCE.' ELSEIF (Cond.EQ.3) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE F(P)' WRITE(9,*)' ' WRITE(9,*)'IS WITHIN THE DESIRED TOLERANCE.' ELSEIF (Cond.EQ.4) THEN WRITE(9,*)'THE APPROXIMATION P AND THE FUNCTION VALUE ' WRITE(9,*)' ' WRITE(9,*)'F(P) ARE BOTH WITHIN THE DESIRED TOLERANCES.' ENDIF WRITE(9,*)' ' RETURN END PROGRAM SECMET C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.6 (Secant Method). C Section 2.4, Newton-Raphson and Secant Methods, Page 85 C PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=100) INTEGER Cond,K REAL Dp,P0,P1,P2,Q0,Q1,RelErr CHARACTER ANS*1 EXTERNAL F 10 CONTINUE CALL INPUTS(P0,P1) Q0=P0 Q1=P1 CALL SECANT(F,Q0,Q1,Delta,Epsilon,Max,P2,Dp,Cond,K) CALL RESULT(P0,P1,P2,Dp,Cond,K) WRITE(9,*)' ' WRITE(9,*) 'WANT TO TRY A ANOTHER STARTING VALUE? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y'.OR.ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) REAL X F=X*X*X-3*X+2 RETURN END SUBROUTINE PRINTFN WRITE(9,*)'F(X) = X*X*X-3*X+2' RETURN END SUBROUTINE SECANT(F,P0,P1,Delta,Epsilon,Max,P2,Dp,Cond,K) PARAMETER(Small=1E-10) INTEGER Cond,K REAL Delta,Epsilon,Df,Dp,P0,P1,P2,Y0,Y1,Y2,RelErr EXTERNAL F K=0 Cond=0 Y0=F(P0) Y1=F(P1) WHILE ((K.LT.Max).AND.(Cond.EQ.0)) Df=(Y1-Y0)/(P1-P0) IF (Df.EQ.0) THEN Cond=1 Dp=P1-P0 P2=P1 ELSE Dp=Y1/Df P2=P1-Dp ENDIF Y2=F(P2) RelErr=ABS(Dp)/(ABS(P2)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y2).LT.Epsilon) Cond=3 IF (Cond.EQ.2).AND.(ABS(Y2).LT.Epsilon) Cond=4 P0=P1 P1=P2 Y0=Y1 Y1=Y2 K=K+1 WRITE(9,1000) K,P1,Y1 REPEAT PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE XSECANT(F,P0,P1,Delta,Epsilon,Max,P2,Dp,Cond,K) C This subroutine uses simulated WHILE loop(s). PARAMETER(Small=1E-10) INTEGER Cond,K REAL Delta,Epsilon,Df,Dp,P0,P1,P2,Y0,Y1,Y2,RelErr EXTERNAL F K=0 Cond=0 Y0=F(P0) Y1=F(P1) 10 IF ((K.LT.Max).AND.(Cond.EQ.0)) THEN Df=(Y1-Y0)/(P1-P0) IF (Df.EQ.0) THEN Cond=1 Dp=P1-P0 P2=P1 ELSE Dp=Y1/Df P2=P1-Dp ENDIF Y2=F(P2) RelErr=ABS(Dp)/(ABS(P2)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y2).LT.Epsilon) Cond=3 IF ((RelErr.LE.Delta).AND.(ABS(Y2).LT.Epsilon)) Cond=4 P0=P1 P1=P2 Y0=Y1 Y1=Y2 K=K+1 WRITE(9,1000) K,P1,Y1 GOTO 10 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE INPUTS(P0,P1) INTEGER I REAL P0,P1 DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE SECANT METHOD IS USED TO FIND A ZERO OF THE FUNCTIO +N:' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'TWO INITIAL APPROXIMATIONS P0, P1 ARE NEEDED.' WRITE(9,*)' ' WRITE(9,*)'ENTER P0 = ' READ(9,*) P0 WRITE(9,*)'ENTER P1 = ' READ(9,*) P1 WRITE(9,*)' ' RETURN END SUBROUTINE RESULT(P0,P1,P2,Dp,Cond,K) INTEGER Cond,I,K REAL Dp,P0,P1,P2,Y2 DO 10 I=1,14 WRITE(9,*)' ' 10 CONTINUE Y2=F(P2) WRITE(9,*)'THE SECANT METHOD WAS USED TO FIND A ZERO OF THE FUNCTI +ON' WRITE(9,*)' ' CALL PRINTFN WRITE(9,*)' ' WRITE(9,*)'STARTING WITH THE TWO APPROXIMATIONS:' WRITE(9,*)' ' WRITE(9,*)'P0 =',P0,' AND P1 =',P1 WRITE(9,*)' ' WRITE(9,*)'AFTER ',K,' ITERATIONS AN APPROXIMATE VALUE FOR THE ZER +O IS:' WRITE(9,*)' ' WRITE(9,*)' P =',P2 WRITE(9,*)' ' WRITE(9,*)' DP =',ABS(Dp),' IS THE ESTIMATED ACCURACY FOR P.' WRITE(9,*)' ' WRITE(9,*)' F(',P2,' ) = ',Y2 WRITE(9,*)' ' IF (Y2.EQ.0) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE IS EXACTLY ZERO!' WRITE(9,*)' ' ENDIF IF (Cond.EQ.0) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE' WRITE(9,*)' ' WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ELSEIF (Cond.EQ.1) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE ' WRITE(9,*)' ' WRITE(9,*)'DIVISION BY ZERO WAS ENCOUNTERED.' ELSEIF (Cond.EQ.2) THEN WRITE(9,*)'THE APPROXIMATION P IS' WRITE(9,*)' ' WRITE(9,*)'WITHIN THE DESIRED TOLERANCE.' ELSEIF (Cond.EQ.3) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE F(P)' WRITE(9,*)' ' WRITE(9,*)'IS WITHIN THE DESIRED TOLERANCE.' ELSEIF (Cond.EQ.4) THEN WRITE(9,*)'THE APPROXIMATION P AND THE FUNCTION VALUE ' WRITE(9,*)' ' WRITE(9,*)'F(P) ARE BOTH WITHIN THE DESIRED TOLERANCES.' ENDIF RETURN END PROGRAM STEFACC C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.7 (Steffensen's Acceleration). C Section 2.5, Aitken's Process & Steffensen's & Muller's Methods, Page 96 C PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=99) INTEGER Cond,K REAL Dp,P0,P3 CHARACTER ANS*1 EXTERNAL F,F1 10 CALL INPUT(P0) CALL STEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K) CALL RESULTS(P0,P3,Dp,Cond,K) WRITE(9,*)' ' WRITE(9,*)'WANT TO TRY ANOTHER STARTING VALUE ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) REAL X F=X*X*X-3*X+2 RETURN END REAL FUNCTION F1(X) REAL X F1=3*X*X-3 RETURN END SUBROUTINE PRINTFUN WRITE(9,*)'F(X) = X*X*X-3*X+2' RETURN END SUBROUTINE STEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K) PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL D1,D2,Delta,DF0,DF1,Dp,Epsilon,P0,P1,P2,P3,Pold,RelErr,Y3 EXTERNAL F,F1 Pold=P0 K=0 Cond=0 P3=P0 P2=P0+1 P1=P0+2 WHILE ((K.LT.Max).AND.(Cond.EQ.0)) P0=P3 DF0=F1(P0) IF (DF0.NE.0) THEN P1=P0 - F(P0)/DF0 ELSE Cond=1 Dp=P3-P2 P3=P0 ENDIF DF1=F1(P1) IF (DF1.EQ.0) THEN Cond=1 Dp=P1-P0 P3=P1 ELSE P2=P1 - F(P1)/DF1 D1=(P1-P0)*(P1-P0) D2=P2-2*P1+P0 IF (D2.EQ.0) THEN Cond=1 Dp=P2-P1 P3=P2 ELSE P3=P0-D1/D2 Dp=P3-P2 ENDIF Y3=F(P3) RelErr=ABS(Dp)/(ABS(P3)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y3).LT.Epsilon) Cond=3 IF (Cond.EQ.2).AND.(ABS(Y3).LT.Epsilon) Cond=4 ENDIF K=K+1 WRITE(9,1000) K,P3,P3 REPEAT P0=Pold PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE XSTEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K) C This subroutine uses simulated WHILE loop(s). PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL D1,D2,Delta,DF0,DF1,Dp,Epsilon,P0,P1,P2,P3,Pold,RelErr,Y3 EXTERNAL F,F1 Pold=P0 K=0 Cond=0 P3=P0 P2=P0+1 P1=P0+2 10 IF ((K.LT.Max).AND.(Cond.EQ.0)) THEN P0=P3 DF0=F1(P0) IF (DF0.NE.0) THEN P1=P0 - F(P0)/DF0 ELSE Cond=1 Dp=P3-P2 P3=P0 ENDIF DF1=F1(P1) IF (DF1.EQ.0) THEN Cond=1 Dp=P1-P0 P3=P1 ELSE P2=P1 - F(P1)/DF1 D1=(P1-P0)*(P1-P0) D2=P2-2*P1+P0 IF (D2.EQ.0) THEN Cond=1 Dp=P2-P1 P3=P2 ELSE P3=P0-D1/D2 Dp=P3-P2 ENDIF Y3=F(P3) RelErr=ABS(Dp)/(ABS(P3)+Small) IF (RelErr.LT.Delta) Cond=2 IF (ABS(Y3).LT.Epsilon) Cond=3 IF (RelErr.LT.Delta).AND.(ABS(Y3).LT.Epsilon) Cond=4 ENDIF K=K+1 WRITE(9,1000) K,P3,P3 GOTO 10 ENDIF P0=Pold PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE INPUT(P0) INTEGER I REAL P0 DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'STEFFENSEN`S ACCELERATION OF THE NEWTON-RAPHSON' WRITE(9,*)' ' WRITE(9,*)'METHOD IS USED TO FIND A ZERO OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'ONE INITIAL APPROXIMATION P0 IS NEEDED.' WRITE(9,*)' ' WRITE(9,*)'ENTER P0 = ' READ(9,*) P0 WRITE(9,*)' ' RETURN END SUBROUTINE RESULTS(P0,P3,Dp,Cond,K) INTEGER Cond,I,K REAL P0,P3,Dp DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'STEFFENSEN`S ACCELERATION OF THE NEWTON-RAPHSON' WRITE(9,*)' ' WRITE(9,*)'METHOD WAS USED TO FIND A ZERO OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'STARTING WITH THE APPROXIMATION P0 =',P0 WRITE(9,*)' ' WRITE(9,*)'AFTER ',K,' ITERATIONS AN APPROXIMATE VALUE OF THE ZERO + IS' WRITE(9,*)' ' WRITE(9,*)' P =',P3 WRITE(9,*)' ' WRITE(9,*)' DP =',ABS(Dp),' IS AN ESTIMATE OF THE ACCURACY.' WRITE(9,*)' ' WRITE(9,*)' F(',P3,' ) =',F(P3) WRITE(9,*)' ' IF (F(P3).EQ.0) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE IS EXACTLY ZERO!' WRITE(9,*)' ' ENDIF IF (Cond.EQ.0) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE' WRITE(9,*)' ' WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ELSEIF (Cond.EQ.1) THEN WRITE(9,*)'CONVERGENCE IS DOUBTFUL. DIVISION BY ZERO WAS ENCOUNT +ERED.' ELSEIF (Cond.EQ.2) THEN WRITE(9,*)'THE APPROXIMATION P IS WITHIN THE DESIRED TOLERANCE.' ELSEIF (Cond.EQ.3) THEN WRITE(9,*)'THE FUNCTION VALUE F(P) IS WITHIN THE DESIRED TOLERAN +CE.' ELSEIF (Cond.EQ.4) THEN WRITE(9,*)'THE APPROXIMATION P AND THE FUNCTION VALUE ' WRITE(9,*)'F(P) ARE BOTH WITHIN THE DESIRED TOLERANCES.' ENDIF RETURN END PROGRAM MULMETH C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.8 (Muller's Method). C Section 2.5, Aitken's Process & Steffensen's & Muller's Methods, Page 97 C PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=99) INTEGER Cond,K REAL P0,P1,P2,P3,Q0,Q1,Q2,Z CHARACTER ANS*1 EXTERNAL F 10 CALL INPUTS(P0,P1,P2) Q0=P0 Q1=P1 Q2=P2 CALL MULLER(F,Q0,Q1,Q2,Delta,Epsilon,Max,P3,Z,K,Cond) CALL RESULT(P0,P1,P2,P3,Z,K,Cond) WRITE(9,*) WRITE(9,*) 'WANT TO TRY NEW STARTING VALUES ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) REAL X F=X*X*X-3*X+2 RETURN END SUBROUTINE PRINTFUN WRITE(9,*)'F(X) = X*X*X-3*X+2' RETURN END SUBROUTINE MULLER(F,P0,P1,P2,Delta,Epsilon,Max,P3,Z,K,Cond) PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL A,B,C,Det,Delta,Disc,Epsilon,H0,H1,E0,E1 REAL P0,P1,P2,P3,Y0,Y1,Y2,RelErr,U,V,Z EXTERNAL F K=0 Cond=0 Y0=F(P0); Y1=F(P1); Y2=F(P2) WHILE ((K.LT.Max).AND.(Cond.EQ.0)) H0=P0-P2 H1=P1-P2 C=Y2 E0=Y0-C E1=Y1-C Det=H0*H1*(H0-H1) A=(E0*H1-H0*E1)/Det B=(H0*H0*E1-H1*H1*E0)/Det IF ((B*B).GT.(4*A*C)) THEN Disc=SQRT(B*B-4*A*C) ELSE DISC=0 ENDIF IF (B.LT.0) Disc=-Disc Z=-2*C/(B+Disc) P3=P2+Z IF (ABS(P3-P1).LT.ABS(P3-P0)) THEN U=P1; P1=P0; P0=U V=Y1; Y1=Y0; Y0=V ENDIF IF (ABS(P3-P2).LT.ABS(P3-P1)) THEN U=P2; P2=P1; P1=U V=Y2; Y2=Y1; Y1=V ENDIF P2=P3; Y2=F(P2) RelErr=ABS(Z)/(ABS(P2)+Small) IF (RelErr.LT.Delta).AND.(ABS(Y2).LT.Epsilon) Cond=1 K=K+1 WRITE(9,1000) K,P2,Y2 REPEAT PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE XMULLER(F,P0,P1,P2,Delta,Epsilon,Max,P3,Z,K,Cond) C This subroutine uses simulated WHILE loop(s). PARAMETER(Small=1E-20) INTEGER Cond,K,Max REAL A,B,C,Det,Delta,Disc,Epsilon,H0,H1,E0,E1 REAL P0,P1,P2,P3,Y0,Y1,Y2,RelErr,U,V,Z EXTERNAL F K=0 Cond=0 Y0=F(P0) Y1=F(P1) Y2=F(P2) 10 IF ((K.LT.Max).AND.(Cond.EQ.0)) THEN H0=P0-P2 H1=P1-P2 C=Y2 E0=Y0-C E1=Y1-C Det=H0*H1*(H0-H1) A=(E0*H1-H0*E1)/Det B=(H0*H0*E1-H1*H1*E0)/Det IF ((B*B).GT.(4*A*C)) THEN Disc=SQRT(B*B-4*A*C) ELSE DISC=0 ENDIF IF (B.LT.0) Disc=-Disc Z=-2*C/(B+Disc) P3=P2+Z IF (ABS(P3-P1).LT.ABS(P3-P0)) THEN U=P1; P1=P0; P0=U V=Y1; Y1=Y0; Y0=V ENDIF IF (ABS(P3-P2).LT.ABS(P3-P1)) THEN U=P2; P2=P1; P1=U V=Y2; Y2=Y1; Y1=V ENDIF P2=P3 Y2=F(P2) RelErr=ABS(Z)/(ABS(P2)+Small) IF (RelErr.LT.Delta).AND.(ABS(Y2).LT.Epsilon) Cond=1 K=K+1 WRITE(9,1000) K,P2,Y2 GOTO 10 ENDIF PAUSE RETURN 1000 FORMAT(I2,4X,F15.7,4X,F15.7) END SUBROUTINE INPUTS(P0,P1,P2) INTEGER I REAL P0,P1,P2 DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'MULLER`S METHOD IS USED TO FIND A ZERO OF THE FUNCTION' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'THREE INITIAL APPROXIMATIONS P0, P1 AND P2 ARE NEEDED.' WRITE(9,*)' ' WRITE(9,*)'ENTER P0 = ' READ(9,*) P0 WRITE(9,*)'ENTER P1 = ' READ(9,*) P1 WRITE(9,*)'ENTER P2 = ' READ(9,*) P2 WRITE(9,*)' ' RETURN END SUBROUTINE RESULT(P0,P1,P2,P3,Z,K,Cond) INTEGER Cond,I,K REAL P0,P1,P2,P3,Y3,Z Y3=F(P3) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'MULLER`S METHOD WAS USED TO FIND A ZERO OF THE FUNCTION +' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'THE THREE INITIAL APPROXIMATIONS WERE:' WRITE(9,*)' ' WRITE(9,*)'P0 =',P0,' P1 =',P1,' P2 =',P2 WRITE(9,*)' ' WRITE(9,*)'AFTER ',K,' ITERATIONS AN APPROXIMATE VALUE FOR A ZERO +IS:' WRITE(9,*)' ' WRITE(9,*)' P =',P3 WRITE(9,*)' ' WRITE(9,*)' F(',P3,' ) = ',Y3 WRITE(9,*)' ' WRITE(9,*)' ',ABS(Z),' IS THE ESTIMATED ACCURACY' WRITE(9,*)' ' IF (Y3.EQ.0) THEN WRITE(9,*)'THE COMPUTED FUNCTION VALUE IS EXACTLY ZERO!' WRITE(9,*)' ' ENDIF IF (Cond.EQ.1) THEN WRITE(9,*)'THE ZERO WAS FOUND AND IS' WRITE(9,*)'WITHIN THE DESIRED TOLERANCES.' ELSE WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE THE' WRITE(9,*)'MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.' ENDIF RETURN END PROGRAM NLSEIDEL C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.9 (Nonlinear Seidel Iteration). C Section 2.6, Iteration for Nonlinear Systems, Page 108 C PARAMETER(Max=199,Tol=1E-6) INTEGER Count REAL P0,Q0,Pold,Qold,Pnew,Qnew,Sep CHARACTER Ans*1 EXTERNAL G1,G2 10 CALL INPUTS(P0,Q0,Pold,Qold) CALL ITERATE(G1,G2,Pold,Qold,Pnew,Qnew,Tol,Max,Sep,Count) CALL RESULTS(P0,Q0,Pold,Qold,Pnew,Qnew,Tol,Sep,Count) WRITE(9,*)' ' WRITE(9,*)'Want to try a different starting point ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 END REAL FUNCTION G1(X,Y) REAL X,Y G1=(X*X - Y + 0.5)/2 RETURN END REAL FUNCTION G2(X,Y) REAL X,Y G2=(- X*X - 4*Y*Y + 8*Y + 4)/8 RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' ' WRITE(9,*)' G (X,Y) = (X*X - Y + 0.5)/2' WRITE(9,*)' 1' WRITE(9,*)' ' WRITE(9,*)' G (X,Y) = (- X*X - 4*Y*Y + 8*Y + 4)/8' WRITE(9,*)' 2' WRITE(9,*)' ' RETURN END SUBROUTINE ITERATE(G1,G2,Pold,Qold,Pnew,Qnew,Tol,Max,Sep,Count) PARAMETER(Big=1E9) INTEGER Count,K,Max REAL P0,Pnew,Pold,Q0,Qnew,Qold,Sep,Tol EXTERNAL G1,G2 Sep=1 K=0 Pnew=Pold Qnew=Qold WHILE (K.LT.Max).AND.(Sep.GT.Tol).AND. & (ABS(Pnew).LT.Big).AND.(ABS(Qnew).LT.Big) Pold=Pnew Qold=Qnew K=K+1 Pnew=G1(Pold,Qold) Qnew=G2(Pnew,Qold) Sep=ABS(Pnew-Pold)+ABS(Qnew-Qold) WRITE(9,1000) K,Pnew,Qnew REPEAT Count=K 1000 FORMAT(I2,4X,5F18.7) PAUSE RETURN END SUBROUTINE XITERATE(G1,G2,Pold,Qold,Pnew,Qnew,Tol,Max,Sep,Count) C This subroutine uses simulated WHILE loop(s). PARAMETER(Big=1E10) INTEGER Count,K,Max REAL P0,Pnew,Pold,Q0,Qnew,Qold,Sep,Tol EXTERNAL G1,G2 Sep=1 K=0 Pnew=Pold Qnew=Qold WRITE(9,999) K,Pnew,Qold 10 IF (K.LT.Max.AND.Sep.GT.Tol) THEN Pold=Pnew Qold=Qnew K=K+1 Pnew=G1(Pold,Qold) Qnew=G2(Pnew,Qold) Sep=ABS(Pnew-Pold)+ABS(Qnew-Qold) WRITE(9,999) K,Pnew,Qold IF (ABS(Pnew).GT.Big).OR.(ABS(Qnew).GT.Big) GOTO 20 GOTO 10 ENDIF 20 CONTINUE Count=K 999 FORMAT(I2,4X,5F18.7) PAUSE RETURN END SUBROUTINE INPUTS(P0,Q0,Pold,Qold) INTEGER I REAL P0,Q0,Pold,Qold DO 10 I=1,10 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' Seidel iteration is used to find the fixed' WRITE(9,*)' ' WRITE(9,*)' point (P,Q) of the non-linear system' WRITE(9,*)' ' WRITE(9,*)' P = G (P,Q), Q = G (P,Q). The functions are:' WRITE(9,*)' 1 2' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ENTER the initial starting point (P ,Q )' WRITE(9,*)' 0 0 ' WRITE(9,*)' ' WRITE(9,*)' P(0) = ' READ(9,*) Pold P0=Pold WRITE(9,*)' ' WRITE(9,*)' Q(0) = ' READ(9,*) Qold Q0=Qold RETURN END SUBROUTINE RESULTS(P0,Q0,Pold,Qold,Pnew,Qnew,Tol,Sep,Count) INTEGER Count,K REAL P0,Pnew,Pold,Q0,Qnew,Qold,Sep,Tol K=Count WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Seidel iteration was used to find the fixed point of th +e' WRITE(9,*)' ' WRITE(9,*)'non-linear system P = G (P,Q), Q = G (P,Q) where:' WRITE(9,*)' 1 2' CALL PRINTFUN WRITE(9,*)'starting with P =',P0,' and Q =',Q0 WRITE(9,*)' 0 0' WRITE(9,*)' ' IF (Sep.LT.Tol) THEN WRITE(9,*)'after ',Count,' iterations Seidel`s method was' WRITE(9,*)' ' WRITE(9,*)'successful and found a solution point (P,Q).' WRITE(9,*)' ' WRITE(9,*)' P =',Pnew,' Q =',Qnew WRITE(9,*)' ' WRITE(9,*)' DP =',Pnew-Pold,' DQ =',Qnew-Qol +d ELSE WRITE(9,*)'Seidel iteration did NOT converge.' WRITE(9,*)' ' WRITE(9,*)'The status after ',Count,' iterations is:' WRITE(9,*)' ' WRITE(9,*)' P(',K,') =',Pnew,' Q(',K,') =',Qnew WRITE(9,*)' ' WRITE(9,*)' DP =',Pnew-Pold,' DQ =',Qnew-Qol +d ENDIF RETURN END PROGRAM NEW2DIM C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 2.10 (Newton-Raphson Method in 2-Dimensions). C Section 2.7, Newton's Method for Systems, Page 116 C PARAMETER(Max=99,Delta=1E-6,Epsilon=1E-6) INTEGER Cond,K REAL D,Pstart,P0,P1,F1 DIMENSION F1(1:2),Pstart(1:2),P0(1:2),P1(1:2) DIMENSION D(1:2,1:2) CHARACTER Ans*1 10 CALL INPUTS(Pstart,P0) CALL NEWTON(P0,P1,F1,D,Cond,K,Delta,Epsilon,Max) CALL RESULTS(Pstart,P0,P1,F1,Cond,K) WRITE(9,*)' ' WRITE(9,*)'Want to try a different starting point ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 END SUBROUTINE FUN(P,F) REAL F,P,X,Y DIMENSION F(1:2),P(1:2) X=P(1) Y=P(2) F(1)=X*X-2*X-Y+0.5 F(2)=X*X+4*Y*Y-4 RETURN END SUBROUTINE JACOBIAN(P,D) REAL D,P,X,Y DIMENSION D(1:2,1:2),P(1:2) X=P(1) Y=P(2) D(1,1)=2*X-2 D(1,2)=-1 D(2,1)=2*X D(2,2)=8*Y RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' ' WRITE(9,*)' 0 = F (X,Y) = X*X - 2*X - Y + 0.5' WRITE(9,*)' 1' WRITE(9,*)' ' WRITE(9,*)' 0 = F (X,Y) = X*X + 4*Y*Y - 4' WRITE(9,*)' 2' WRITE(9,*)' ' RETURN END SUBROUTINE NEWTON(P0,P1,F1,D,Cond,K,Delta,Epsilon,Max) INTEGER Cond,K,Max REAL D,Delta,Det,DP,Epsilon,F0,F1,FnZero,P0,P1,Small,RelErr DIMENSION DP(1:2),F0(1:2),F1(1:2),Pstart(1:2) DIMENSION P0(1:2),P1(1:2),D(1:2,1:2) Small=Delta K=0 Cond=0 CALL FUN(P0,F0) P1(1)=P0(1) P1(2)=P0(2) F1(1)=F0(1) F1(2)=F0(2) WRITE(9,1000) K,( P1(J), J=1,2 ) WHILE (K.LT.Max.AND.Cond.EQ.0) P0(1)=P1(1) P0(2)=P1(2) F0(1)=F1(1) F0(2)=F1(2) K=K+1 CALL JACOBIAN(P0,D) Det=D(1,1)*D(2,2)-D(1,2)*D(2,1) IF (Det.EQ.0) THEN DP(1)=0 DP(2)=0 Cond=1 ELSE DP(1)=(F0(1)*D(2,2)-F0(2)*D(1,2))/Det DP(2)=(F0(2)*D(1,1)-F0(1)*D(2,1))/Det ENDIF P1(1)=P0(1)-DP(1) P1(2)=P0(2)-DP(2) CALL FUN(P1,F1) RelErr=(ABS(DP(1))+ABS(DP(2)))/(ABS(P1(1))+ABS(P1(2))+Small) FnZero=ABS(F1(1))+ABS(F1(2)) IF (RelErr.LT.Delta.AND.FnZero.LT.Epsilon) THEN IF (Cond.NE.1) Cond=2 ENDIF WRITE(9,1000) K, ( P1(J), J=1,2 ) REPEAT 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE XNEWTON(P0,P1,F1,D,Cond,K,Delta,Epsilon,Max) C This subroutine uses simulated WHILE loop(s). INTEGER Cond,K,Max REAL D,Delta,Det,DP,Epsilon,F0,F1,FnZero,P0,P1,Small,RelErr DIMENSION DP(1:2),F0(1:2),F1(1:2),Pstart(1:2) DIMENSION P0(1:2),P1(1:2),D(1:2,1:2) Small=Delta K=0 Cond=0 CALL FUN(P0,F0) P1(1)=P0(1) P1(2)=P0(2) F1(1)=F0(1) F1(2)=F0(2) WRITE(9,1000) K,( P1(J), J=1,2 ) 10 IF (K.LT.Max.AND.Cond.EQ.0) THEN P0(1)=P1(1) P0(2)=P1(2) F0(1)=F1(1) F0(2)=F1(2) K=K+1 CALL JACOBIAN(P0,D) Det=D(1,1)*D(2,2)-D(1,2)*D(2,1) IF (Det.EQ.0) THEN DP(1)=0 DP(2)=0 Cond=1 ELSE DP(1)=(F0(1)*D(2,2)-F0(2)*D(1,2))/Det DP(2)=(F0(2)*D(1,1)-F0(1)*D(2,1))/Det ENDIF P1(1)=P0(1)-DP(1) P1(2)=P0(2)-DP(2) CALL FUN(P1,F1) RelErr=(ABS(DP(1))+ABS(DP(2)))/(ABS(P1(1))+ABS(P1(2))+Small) FnZero=ABS(F1(1))+ABS(F1(2)) IF (RelErr.LT.Delta.AND.FnZero.LT.Epsilon) THEN IF (Cond.NE.1) Cond=2 ENDIF WRITE(9,1000) K, ( P1(J), J=1,2 ) GOTO 10 ENDIF 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE INPUTS(Pstart,P0) REAL Pstart,P0 DIMENSION Pstart(1:2),P0(1:2) DO 10 I=1,10 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' Newton-Raphson iteration is used to find a' WRITE(9,*)' ' WRITE(9,*)' solution point (P,Q) of the non-linear system' WRITE(9,*)' ' WRITE(9,*)' 0 = F (X,Y)' WRITE(9,*)' 1' WRITE(9,*)' ' WRITE(9,*)' 0 = F (X,Y) The current system of equations is:' WRITE(9,*)' 2' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ENTER the initial starting point (P ,Q )' WRITE(9,*)' 0 0 ' WRITE(9,*)' ' WRITE(9,*)' P(0) = ' READ(9,*) P0(1) Pstart(1)=P0(1) WRITE(9,*)' ' WRITE(9,*)' Q(0) = ' READ(9,*) P0(2) Pstart(2)=P0(2) RETURN END SUBROUTINE RESULTS(Pstart,P0,P1,F1,Cond,K) INTEGER Cond,K REAL F1,P0,P1 DIMENSION F1(1:2),Pstart(1:2),P0(1:2),P1(1:2) K=Count WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Newton-Raphson iteration was used to solve the non-line +ar system' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)'starting with P =',Pstart(1),' and Q =',Pstart(2) WRITE(9,*)' 0 0' WRITE(9,*)' ' WRITE(9,*)'It took ',K,' iterations to compute the solution point +(P,Q).' WRITE(9,*)' ' IF (Cond.EQ.0) THEN WRITE(9,*)'However, the maximum number of iterations was exceede +d.' ELSEIF (Cond.EQ.1) THEN WRITE(9,*)'However, division by zero was encountered.' ELSEIF (Cond.EQ.2) THEN WRITE(9,*)'The solution is within the desired tolerances.' ENDIF WRITE(9,*)' ' IF (Cond.EQ.0.OR.Cond.EQ.1) THEN WRITE(9,*)' P(',K,') =',P1(1),' , Q(',K,') =',P1(2) ELSEIF (Cond.EQ.2) THEN WRITE(9,*)' P =',P1(1),' , Q =',P1(2) ENDIF WRITE(9,*)' ' WRITE(9,*)' |DP| =',ABS(P1(1)-P0(1)),' , |DQ| =',ABS(P +1(2)-P0(2)) WRITE(9,*)' ' WRITE(9,*)' F (P,Q) =',F1(1) WRITE(9,*)' 1' WRITE(9,*)' F (P,Q) =',F1(2) WRITE(9,*)' 2' RETURN END