PROGRAM SOLVEDE 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 9.1 (Euler's Method). C Section 9.2, Euler's Method, Page 435 C C Algorithm 9.2 (Heun's Method). C Section 9.3, Heun's Method, Page 441 C C Algorithm 9.3 (Taylor's Method of Order 4). C Section 9.4, Taylor Series Method, Page 448 C Algorithm 9.4 (Runge-Kutta Method of Order 4). C Section 9.5, Runge-Kutta Methods, Page 460 C PARAMETER(MaxM=1000) INTEGER DoMo,M,Mend,Meth,Order,State REAL A,B,D,T,TK,Y,YK,Y0 DIMENSION D(1:4),T(0:1000),Y(0:1000) CHARACTER Ans*1 EXTERNAL DF,F Meth=1 A=0 B=0 Y0=0 M=1 State=1 10 CONTINUE CALL MESSAGE(Meth,Order,State) DoMo=1 20 CONTINUE CALL INPUT(Meth) 30 CONTINUE CALL EPOINTS(A,B,Y0,M,State) IF (Meth.EQ.1) THEN CALL EULER(F,A,B,Y0,M,Mend,T,Y) ELSEIF (Meth.EQ.2) THEN CALL CAUCHY(F,A,B,Y0,M,Mend,T,Y) ELSEIF (Meth.EQ.3) THEN CALL HEUN(F,A,B,Y0,M,Mend,T,Y) ELSEIF (Meth.EQ.4) THEN CALL TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order) ELSEIF (Meth.EQ.5) THEN CALL TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order) ELSEIF (Meth.EQ.6) THEN CALL RUNGE(F,A,B,Y0,M,Mend,T,Y) ENDIF CALL RESULTS(T,Y,M,Mend,Meth) WRITE(9,*)' ' WRITE(9,*)' Want to try another initial condition ? & ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN State=2 ELSE State=0 ENDIF IF (State.EQ.0.OR.State.EQ.1) GOTO 30 WRITE(9,*)' Want to change the differential equation ? & ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN DoMo=0 ELSE State=0 ENDIF IF (DoMo.EQ.1) GOTO 20 WRITE(9,*)'Want to try another method of approximation ? ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN Meth=0 ELSE State=0 ENDIF IF (Meth.NE.0) GOTO 10 RETURN END REAL FUNCTION F(T,Y) REAL T,Y F=T**2+Y**2 RETURN END SUBROUTINE DF(T,Y,D) REAL D,T,Y DIMENSION D(1:4) D(1) = ( T - Y )/2 D(2) = ( 2 - T + Y)/4 D(3) = (-2 + T - Y)/8 D(4) = ( 2 - T + Y)/16 RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' Y` = T**2+Y**2' RETURN END SUBROUTINE EULER(F,A,B,Y0,M,Mend,T,Y) PARAMETER(Big=1E13) INTEGER K,M,Mend REAL A,B,H,T,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO K=0,M-1 Y(K+1)=Y(K)+H*F(T(K),Y(K)) T(K+1)=A+H*(K+1) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN ENDDO RETURN END SUBROUTINE HEUN(F,A,B,Y0,M,Mend,T,Y) PARAMETER(Big=1E13) INTEGER J,M,Mend REAL A,B,H,K1,K2,P,T,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO J=0,M-1 K1=F(T(J),Y(J)) P=Y(J)+H*K1 T(J+1)=A+H*(J+1) Mend=J+1 IF (Big.LT.ABS(P)) THEN Y(J+1)=P RETURN ENDIF K2=F(T(J+1),P) Y(J+1)=Y(J)+H*(K1+K2)/2 IF (Big.LT.ABS(Y(J+1))) RETURN ENDDO RETURN END SUBROUTINE CAUCHY(F,A,B,Y0,M,Mend,T,Y) PARAMETER(Big=1E13) INTEGER J,M,Mend REAL A,B,H,K1,K2,P,T,TP,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO J=0,M-1 K1=F(T(J),Y(J)) P=Y(J)+H*K1/2 T(J+1)=A+H*(J+1) Mend=J+1 IF (Big.LT.ABS(P)) THEN Y(J+1)=P RETURN ENDIF TP=T(J)+H/2 K2=F(TP,P) Y(J+1)=Y(J)+H*K2 IF (Big.LT.ABS(Y(J+1))) RETURN ENDDO RETURN END SUBROUTINE TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order) PARAMETER(Big=1E13) INTEGER J,M,Mend,Order REAL A,B,D,H,T,TK,Y,YK,Y0 DIMENSION D(1:4),T(0:1000),Y(0:1000) EXTERNAL DF H=(B-A)/M T(0)=A Y(0)=Y0 DO K=0,M-1 TK=T(K) YK=Y(K) CALL DF(TK,YK,D) IF (Order.EQ.3) THEN Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*D(3)/6)) ELSE Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*(D(3)/6+H*D(4)/24))) ENDIF T(K+1)=A+H*(K+1) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN ENDDO RETURN END SUBROUTINE RUNGE(F,A,B,Y0,M,Mend,T,Y) PARAMETER(Big=1E13) INTEGER J,M,Mend,Pole REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO J=0,M-1 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) IF (Big.LT.ABS(YJ+0.5*K1)) RETURN K2=H*F(TJ+H/2,YJ+0.5*K1) IF (Big.LT.ABS(YJ+0.5*K2)) RETURN K3=H*F(TJ+H/2,YJ+0.5*K2) IF (Big.LT.ABS(YJ+K3)) RETURN K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) Mend=J+1 IF (Big.LT.ABS(Y(J+1))) RETURN ENDDO RETURN END SUBROUTINE XRUNGE(F,A,B,Y0,M,Mend,T,Y) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,M,Mend,Pole REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F Pole=0 H=(B-A)/M T(0)=A Y(0)=Y0 DO J=0,M-1 TJ=T(J) YJ=Y(J) Mend=J+1 K1=H*F(TJ,YJ) IF (Big.LT.ABS(YJ+0.5*K1)) THEN Pole=1 GOTO 10 ENDIF K2=H*F(TJ+H/2,YJ+0.5*K1) IF (Big.LT.ABS(YJ+0.5*K2)) THEN Pole=1 GOTO 10 ENDIF K3=H*F(TJ+H/2,YJ+0.5*K2) IF (Big.LT.ABS(YJ+K3)) THEN Pole=1 GOTO 10 ENDIF K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) IF (Big.LT.ABS(Y(J+1))) GOTO 10 ENDDO 10 CONTINUE IF (Pole.EQ.0) RETURN P=YJ+K1 T(J+1)=A+H*(J+1) IF (Big.LT.ABS(P)) THEN Y(J+1)=P RETURN ENDIF K2=F(T(J+1),P) Y(J+1)=YJ+H*(K1+K2)/2 RETURN END SUBROUTINE XEULER(F,A,B,Y0,M,Mend,T,Y) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER K,M,Mend REAL A,B,H,T,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO 10 K=0,M-1 Y(K+1)=Y(K)+H*F(T(K),Y(K)) T(K+1)=A+H*(K+1) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) GOTO 20 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XHEUN(F,A,B,Y0,M,Mend,T,Y) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,M,Mend REAL A,B,H,K1,K2,P,T,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO 10 J=0,M-1 K1=F(T(J),Y(J)) P=Y(J)+H*K1 T(J+1)=A+H*(J+1) Mend=J+1 IF (Big.LT.ABS(P)) THEN Y(J+1)=P GOTO 20 ENDIF K2=F(T(J+1),P) Y(J+1)=Y(J)+H*(K1+K2)/2 IF (Big.LT.ABS(Y(J+1))) GOTO 20 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XCAUCHY(F,A,B,Y0,M,Mend,T,Y) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,M,Mend REAL A,B,H,K1,K2,P,T,TP,Y,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 DO 10 J=0,M-1 K1=F(T(J),Y(J)) P=Y(J)+H*K1/2 T(J+1)=A+H*(J+1) Mend=J+1 IF (Big.LT.ABS(P)) THEN Y(J+1)=P GOTO 20 ENDIF TP=T(J)+H/2 K2=F(TP,P) Y(J+1)=Y(J)+H*K2 IF (Big.LT.ABS(Y(J+1))) GOTO 20 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XTAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,M,Mend,Order REAL A,B,D,H,T,TK,Y,YK,Y0 DIMENSION D(1:4),T(0:1000),Y(0:1000) EXTERNAL DF H=(B-A)/M T(0)=A Y(0)=Y0 DO 10 K=0,M-1 TK=T(K) YK=Y(K) CALL DF(TK,YK,D) IF (Order.EQ.3) THEN Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*D(3)/6)) ELSE Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*(D(3)/6+H*D(4)/24))) ENDIF T(K+1)=A+H*(K+1) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) GOTO 20 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XRUNGE(F,A,B,Y0,M,Mend,T,Y) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,M,Mend,Pole REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0 DIMENSION T(0:1000),Y(0:1000) EXTERNAL F Pole=0 H=(B-A)/M T(0)=A Y(0)=Y0 DO 10 J=0,M-1 TJ=T(J) YJ=Y(J) Mend=J+1 K1=H*F(TJ,YJ) IF (Big.LT.ABS(YJ+0.5*K1)) THEN Pole=1 GOTO 20 ENDIF K2=H*F(TJ+H/2,YJ+0.5*K1) IF (Big.LT.ABS(YJ+0.5*K2)) THEN Pole=1 GOTO 20 ENDIF K3=H*F(TJ+H/2,YJ+0.5*K2) IF (Big.LT.ABS(YJ+K3)) THEN Pole=1 GOTO 20 ENDIF K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) IF (Big.LT.ABS(Y(J+1))) GOTO 20 10 CONTINUE 20 CONTINUE IF (Pole.EQ.0) GOTO 30 P=YJ+K1 T(J+1)=A+H*(J+1) IF (Big.LT.ABS(P)) THEN Y(J+1)=P GOTO 30 ENDIF K2=F(T(J+1),P) Y(J+1)=YJ+H*(K1+K2)/2 30 CONTINUE RETURN END SUBROUTINE MESSAGE(Meth,Order,State) INTEGER K,Meth,Order,State REAL Resp DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' Solution of the differential equation Y` = F +(T,Y)' WRITE(9,*)' ' WRITE(9,*)' with the initial condition Y(A) = Y .' WRITE(9,*)' 0' WRITE(9,*)' ' WRITE(9,*)' A numerical approximation is computed over [A +,B].' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Choose the method of approximation:' WRITE(9,*)' ' WRITE(9,*)' < 1 > Euler`s method' WRITE(9,*)' ' WRITE(9,*)' < 2 > Modified Euler-Cauchy method' WRITE(9,*)' ' WRITE(9,*)' < 3 > Heun`s method' WRITE(9,*)' ' WRITE(9,*)' < 4 > Taylor`s method of order N=3' WRITE(9,*)' ' WRITE(9,*)' < 5 > Taylor`s method of order N=4' WRITE(9,*)' ' WRITE(9,*)' < 6 > Runge-Kutta method of order N=4' WRITE(9,*)' ' WRITE(9,*)' SELECT < 1 - 6 > ? ' Resp=Meth READ(9,*) Resp Meth=INT(Resp) IF ((Meth.LT.1.OR.Meth.GT.6).AND.State.NE.0) THEN Meth=1 ENDIF IF (Meth.EQ.4) THEN Order=3 ENDIF IF (Meth.EQ.5) THEN Order=4 ENDIF RETURN END SUBROUTINE INPUT(Meth) INTEGER Meth REAL Valu CHARACTER Ans*1 DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' You chose ' WRITE(9,*)' ' IF (Meth.EQ.1) THEN WRITE(9,*)' Euler`s method to solve the D.E.' ELSEIF (Meth.EQ.2) THEN WRITE(9,*)' the Modified Euler-Cauchy method to solve + the D.E.' ELSEIF (Meth.EQ.3) THEN WRITE(9,*)' Heun`s method to solve the D.E.' ELSEIF (Meth.EQ.4) THEN WRITE(9,*)' Taylor`s method of order N=3 to solve + the D.E.' ELSEIF (Meth.EQ.5) THEN WRITE(9,*)' Taylor`s method of order N=4 to solve + the D.E.' ELSEIF (Meth.EQ.6) THEN WRITE(9,*)' the Runge-Kutta method of order N=4 to solve + the D.E.' ENDIF WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' with the initial condition Y(A) = Y0.' WRITE(9,*)' ' WRITE(9,*)' A numerical approximation is computed over [A,B].' WRITE(9,*)' ' WRITE(9,*)' You must enter the endpoints for the interval,' WRITE(9,*)' ' WRITE(9,*)' the initial condition Y0, and the number of steps + M.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,'(A)') Ans RETURN END SUBROUTINE EPOINTS(A,B,Y0,M,State) INTEGER I,M,Stat,State REAL A,B,Y0,Valu C STATUS = (Change,Enter,Done) CHARACTER Ans*1 Stat=1 IF (State.EQ.0) THEN Stat=0 ENDIF 10 CONTINUE DO 20 I=1,15 WRITE(9,*)' ' 20 CONTINUE CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' IF (Stat.EQ.1) THEN WRITE(9,*)' ENTER the left endpoint A = ' READ(9,*) A WRITE(9,*)' ' WRITE(9,*)' ENTER the right endpoint B = ' READ(9,*) B WRITE(9,*)' ' WRITE(9,*)' ENTER initial condition Y(A) = ' READ(9,*) Y0 WRITE(9,*)' ' WRITE(9,*)' ENTER the number of steps M = ' Valu=M READ(9,*) Valu M=INT(Valu) IF (M.LT.1) THEN M=1 ENDIF IF (M.GT.1000) THEN M=1000 ENDIF ELSE WRITE(9,*)' The left endpoint is A =',A WRITE(9,*)' ' WRITE(9,*)' The right endpoint is B =',B WRITE(9,*)' ' WRITE(9,*)' Initial condition is Y(A) =',Y0 WRITE(9,*)' ' WRITE(9,*)' The number of steps is M = ',M ENDIF WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Do you want to make a change ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y'.OR.Ans.EQ.'y') THEN Stat=0 DO 30 I=1,15 WRITE(9,*)' ' 30 CONTINUE WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' The current left endpoint is A =',A WRITE(9,*)' ENTER the NEW left endpoint A = ' READ(9,*) A WRITE(9,*)' ' WRITE(9,*)' The current right endpoint is B =',B WRITE(9,*)' ENTER the NEW right endpoint B = ' READ(9,*) B WRITE(9,*)' ' WRITE(9,*)' The current I. C. is Y(A) =',Y0 WRITE(9,*)' ENTER the NEW I. C. Y(A) = ' READ(9,*) Y0 WRITE(9,*)' ' WRITE(9,*)' The current value of M is M = ',M WRITE(9,*)' ENTER the NEW value of M = ' Valu=M READ(9,*) Valu M=INT(Valu) IF (M.LT.1) THEN M=1 ENDIF IF (M.GT.1000) THEN M=1000 ENDIF ELSE Stat=2 ENDIF IF (Stat.EQ.0.OR.Stat.EQ.1) GOTO 10 RETURN END SUBROUTINE RESULTS(T,Y,M,Mend,Meth) INTEGER I,K,M,Mend,Meth REAL T,Y DIMENSION T(0:1000),Y(0:1000) DO 10 I=1,10 WRITE(9,*)' ' 10 CONTINUE IF (Meth.EQ.1) THEN WRITE(9,*)'Euler`s method was used to solve the D.E.' ELSEIF (Meth.EQ.2) THEN WRITE(9,*)'The Modified Euler-Cauchy method was used to solve th +e D.E.' ELSEIF (Meth.EQ.3) THEN WRITE(9,*)'Heun`s method was used to solve the D.E.' ELSEIF (Meth.EQ.4) THEN WRITE(9,*)'Taylor`s method of order N=3 was used to solve the D. +E.' ELSEIF (Meth.EQ.5) THEN WRITE(9,*)'Taylor`s method of order N=4 was used to solve the D. +E.' ELSEIF (Meth.EQ.6) THEN WRITE(9,*)'The Runge-Kutta method of order N=4 was used to solve + the D.E.' ENDIF WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'with Y(',T(0),' ) =',Y(0) WRITE(9,*)' ' WRITE(9,*)' K T(K) Y(K)' WRITE(9,*)' ---------------------------------------' DO 20 K=0,Mend WRITE(9,999) K, T(K), Y(K) 20 CONTINUE IF (Mend.LT.M) THEN WRITE(9,*)' ' WRITE(9,*)'The solution points are approaching a pole.' ENDIF 999 FORMAT(I4,4X,4F15.7) PAUSE RETURN END PROGRAM SOLVEDE 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 9.5 (Runge-Kutta-Fehlberg Method RKF45). C Section 9.5, Runge-Kutta Methods, Page 461 C C Algorithm 9.6 (Adams-Bashforth-Moulton Method). C Section 9.6, Predictor-Corrector Method, Page 471 C C Algorithm 9.7 (Milne-Simpson Method). C Section 9.6, Predictor-Corrector Method, Page 472 C Algorithm 9.8 (The Hamming Method). C Section 9.6, Predictor-Corrector Method, Page 473 C PARAMETER(MaxM=1000) INTEGER DoMo,M,Mend,Meth,Order,State REAL A,B,D,T,TK,Tol,Y,YK,Y0 DIMENSION D(1:4),T(1:1000),Y(1:1000) CHARACTER Ans*1 EXTERNAL DF,F Meth=1 A=0 B=0 Y0=0 M=1 State=1 10 CONTINUE CALL MESSAGE(Meth,Order,State) DoMo=1 20 CONTINUE CALL INPUT(Meth) 30 CONTINUE CALL EPOINTS(A,B,Y0,M,Tol,State) IF (Meth.EQ.1) THEN CALL ABM(F,A,B,Y0,M,T,Y,Mend) ELSEIF (Meth.EQ.2) THEN CALL MILNE(F,A,B,Y0,M,T,Y,Mend) ELSEIF (Meth.EQ.3) THEN CALL HAMMING(F,A,B,Y0,M,T,Y,Mend) ELSEIF (Meth.EQ.6) THEN CALL RKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend) ENDIF CALL RESULTS(T,Y,M,Mend,Meth) WRITE(9,*)' ' WRITE(9,*)' Want to try another initial condition ? ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN State=2 ELSE State=0 ENDIF IF (State.EQ.0.OR.State.EQ.1) GOTO 30 WRITE(9,*)' Want to change the differential equation ? + ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN DoMo=0 ELSE State=0 ENDIF IF (DoMo.EQ.1) GOTO 20 WRITE(9,*)'Want to try another method of approximation ? ' READ(9,'(A)') Ans IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN Meth=0 ELSE State=0 ENDIF IF (Meth.NE.0) GOTO 10 RETURN END REAL FUNCTION F(T,Y) REAL T,Y F=T**2+Y**2 RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' Y` = T**2+Y**2' RETURN END SUBROUTINE HAMMING(F,A,B,Y0,M,T,Y,Mend) PARAMETER(Big=1E13) REAL A,B,Cnew,Cold,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 !Use Runge-Kutta to compute three starting values. DO J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) ENDDO !Start of Procedure Hamming F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) Pold=0 Cold=0 DO K=3,M-1 Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3 Pmod=Pnew+112*(Cold-Pold)/121 T(K+1)=A+H*(K+1) F4=F(T(K+1),Pmod) Cnew=(9*Y(K)-Y(K-2)+3*H*(-F2+2*F3+F4))/8 Y(K+1)=Cnew+9*(Pnew-Cnew)/121 Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN Pold=Pnew Cold=Cnew F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) ENDDO RETURN END SUBROUTINE ABM(F,A,B,Y0,M,T,Y,Mend) PARAMETER(Big=1E13) INTEGER J,K,M REAL A,B,H,H2,F0,F1,F2,F3,F4,P,T,TK,Y,YK,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 !Use Runge-Kutta to compute three starting values. DO J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) ENDDO !Start of Procedure Adams-Bashforth-Moulton F0=F(T(0),Y(0)) F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) H2=H/24 DO K=3,M-1 TK=T(K) YK=Y(K) P=YK+H2*(-9*F0+37*F1-59*F2+55*F3) T(K+1)=A+H*(K+1) F4=F(T(K+1),P) Y(K+1)=YK+H2*(F1-5*F2+19*F3+9*F4) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN F0=F1 F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) ENDDO RETURN END SUBROUTINE MILNE(F,A,B,Y0,M,T,Y,Mend) PARAMETER(Big=1E13) INTEGER J,K,M REAL A,B,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Yold,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 !Use Runge-Kutta to compute three starting values. DO J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) ENDDO !Start of Procedure Milne-Simpson F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) Pold=0 Yold=0 DO K=3,M-1 Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3 Pmod=Pnew+28*(Yold-Pold)/29 T(K+1)=A+H*(K+1) F4=F(T(K+1),Pmod) Y(K+1)=Y(K-1)+H*(F2+4*F3+F4)/3 Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN Pold=Pnew Yold=Y(K+1) F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) ENDDO RETURN END SUBROUTINE RKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend) PARAMETER(Big=1E13) INTEGER I,J,K,M,Mend REAL A,B,T,Tol,Y,Y0 REAL Br,H,Err,Hmin,Hmax,K1,K2,K3,K4,K5,K6,TJ,YJ REAL A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,S REAL A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5 REAL M1,M3,M4,M5,M6,Y1,Y2,Y3,Y4,Y5,Y6,Ygood,Ynew DIMENSION T(1:1000),Y(1:1000) EXTERNAL F DATA A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5, A A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5 B/ 0.25, 0.25, 0.375, 0.09375, 0.28125, 0.9230769231, C 0.8793809741, -3.277196177, 3.320892126, 1.0, D 2.032407407, -8.0, 7.173489279, -0.2058966862, 0.5, E -0.2962962963, 2.0, -1.381676413, 0.4529727096, -0.275, F 0.002777777778, -0.02994152047, -0.02919989367, 0.02, G 0.03636363636, 0.1157407407, 0.5489278752, 0.535331384, -0.2/ H=(B-A)/M Hmin=H/1024 Hmax=H*64 T(0)=A Y(0)=Y0 T(0)=A J=0 TJ=A Br=B-0.00001*ABS(B) WHILE (T(J).LT.B) IF ((T(J)+H).GT.Br) H=B-T(J) TJ=T(J) YJ=Y(J) Y1=YJ K1=H*F(TJ,Y1) Y2=YJ+B2*K1; IF Big.LT.ABS(Y2) RETURN K2=H*F(TJ+A2*H,Y2) Y3=YJ+B3*K1+C3*K2; IF Big.LT.ABS(Y3) RETURN K3=H*F(TJ+A3*H,Y3) Y4=YJ+B4*K1+C4*K2+D4*K3; IF Big.LT.ABS(Y4) RETURN K4=H*F(TJ+A4*H,Y4) Y5=YJ+B5*K1+C5*K2+D5*K3+E5*K4; IF Big.LT.ABS(Y5) RETURN K5=H*F(TJ+A5*H,Y5) Y6=YJ+B6*K1+C6*K2+D6*K3+E6*K4+F6*K5; IF Big.LT.ABS(Y6) RETURN K6=H*F(TJ+A6*H,Y6) Err=ABS(R1*K1+R3*K3+R4*K4+R5*K5+R6*K6) Ynew=YJ+N1*K1+N3*K3+N4*K4+N5*K5 Err=ABS(Err) IF (Err.LT.Tol.OR.H.LE.2*Hmin) THEN Y(J+1)=Ynew IF ((TJ+H).GT.Br) THEN T(J+1)=B ELSE T(J+1)=TJ+H ENDIF J=J+1 TJ=T(J) ENDIF IF (Err.EQ.0) THEN S=0 ELSE S=0.84*(Tol*H/Err)**0.25 ENDIF IF (S.LT.0.75.AND.H.GT.(2*Hmin)) H=H/2 IF (S.GT.1.50.AND.(2*H).LT.Hmax) H=H*2 IF (Big.LT.ABS(Y(J)).OR.MaxM.EQ.J) RETURN Mend=J IF (B.GT.T(J)) THEN M=J+1 ELSE M=J ENDIF WRITE(9,1000) T(J),Y(J) REPEAT RETURN 1000 FORMAT(5X,2F18.9) END SUBROUTINE XHAMMING(F,A,B,Y0,M,T,Y,Mend) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) REAL A,B,Cnew,Cold,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 C Use Runge-Kutta to compute three starting values. DO 10 J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) 10 CONTINUE C Start of Procedure Hamming F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) Pold=0 Cold=0 DO 20 K=3,M-1 Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3 Pmod=Pnew+112*(Cold-Pold)/121 T(K+1)=A+H*(K+1) F4=F(T(K+1),Pmod) Cnew=(9*Y(K)-Y(K-2)+3*H*(-F2+2*F3+F4))/8 Y(K+1)=Cnew+9*(Pnew-Cnew)/121 Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN Pold=Pnew Cold=Cnew F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) 20 CONTINUE RETURN END SUBROUTINE XABM(F,A,B,Y0,M,T,Y,Mend) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,K,M REAL A,B,H,H2,F0,F1,F2,F3,F4,P,T,TK,Y,YK,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 C Use Runge-Kutta to compute three starting values. DO 10 J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) 10 CONTINUE C Start of Procedure Adams-Bashforth-Moulton F0=F(T(0),Y(0)) F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) H2=H/24 DO 20 K=3,M-1 TK=T(K) YK=Y(K) P=YK+H2*(-9*F0+37*F1-59*F2+55*F3) T(K+1)=A+H*(K+1) F4=F(T(K+1),P) Y(K+1)=YK+H2*(F1-5*F2+19*F3+9*F4) Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN F0=F1 F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) 20 CONTINUE RETURN END SUBROUTINE XMILNE(F,A,B,Y0,M,T,Y,Mend) C This subroutine uses labeled DO loop(s). PARAMETER(Big=1E13) INTEGER J,K,M REAL A,B,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Yold,Y0 REAL K1,K2,K3,K4,TJ,YJ DIMENSION T(1:1000),Y(1:1000) EXTERNAL F H=(B-A)/M T(0)=A Y(0)=Y0 C Use Runge-Kutta to compute three starting values. DO 10 J=0,2 TJ=T(J) YJ=Y(J) K1=H*F(TJ,YJ) K2=H*F(TJ+H/2,YJ+0.5*K1) K3=H*F(TJ+H/2,YJ+0.5*K2) K4=H*F(TJ+H,YJ+K3) Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6 T(J+1)=A+H*(J+1) 10 CONTINUE C Start of Procedure Milne-Simpson F1=F(T(1),Y(1)) F2=F(T(2),Y(2)) F3=F(T(3),Y(3)) Pold=0 Yold=0 DO 20 K=3,M-1 Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3 Pmod=Pnew+28*(Yold-Pold)/29 T(K+1)=A+H*(K+1) F4=F(T(K+1),Pmod) Y(K+1)=Y(K-1)+H*(F2+4*F3+F4)/3 Mend=K+1 IF (Big.LT.ABS(Y(K+1))) RETURN Pold=Pnew Yold=Y(K+1) F1=F2 F2=F3 F3=F(T(K+1),Y(K+1)) 20 CONTINUE RETURN END SUBROUTINE XRKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend) PARAMETER(Big=1E13) INTEGER I,J,K,M,Mend REAL A,B,T,Tol,Y,Y0 REAL Br,H,Err,Hmin,Hmax,K1,K2,K3,K4,K5,K6,TJ,YJ REAL A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,S REAL A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5 REAL M1,M3,M4,M5,M6,Y1,Y2,Y3,Y4,Y5,Y6,Ygood,Ynew DIMENSION T(1:1000),Y(1:1000) EXTERNAL F DATA A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5, A A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5 B/ 0.25, 0.25, 0.375, 0.09375, 0.28125, 0.9230769231, C 0.8793809741, -3.277196177, 3.320892126, 1.0, D 2.032407407, -8.0, 7.173489279, -0.2058966862, 0.5, E -0.2962962963, 2.0, -1.381676413, 0.4529727096, -0.275, F 0.002777777778, -0.02994152047, -0.02919989367, 0.02, G 0.03636363636, 0.1157407407, 0.5489278752, 0.535331384, -0.2/ H=(B-A)/M Hmin=H/1024 Hmax=H*64 T(0)=A Y(0)=Y0 T(0)=A J=0 TJ=A Br=B-0.00001*ABS(B) 10 IF (T(J).LT.B) THEN IF ((T(J)+H).GT.Br) H=B-T(J) TJ=T(J) YJ=Y(J) Y1=YJ K1=H*F(TJ,Y1) Y2=YJ+B2*K1 IF (Big.LT.ABS(Y2)) RETURN K2=H*F(TJ+A2*H,Y2) Y3=YJ+B3*K1+C3*K2 IF (Big.LT.ABS(Y3)) RETURN K3=H*F(TJ+A3*H,Y3) Y4=YJ+B4*K1+C4*K2+D4*K3 IF (Big.LT.ABS(Y4)) RETURN K4=H*F(TJ+A4*H,Y4) Y5=YJ+B5*K1+C5*K2+D5*K3+E5*K4 IF (Big.LT.ABS(Y5)) RETURN K5=H*F(TJ+A5*H,Y5) Y6=YJ+B6*K1+C6*K2+D6*K3+E6*K4+F6*K5 IF (Big.LT.ABS(Y6)) RETURN K6=H*F(TJ+A6*H,Y6) Err=ABS(R1*K1+R3*K3+R4*K4+R5*K5+R6*K6) Ynew=YJ+N1*K1+N3*K3+N4*K4+N5*K5 Err=ABS(Err) IF (Err.LT.Tol.OR.H.LE.2*Hmin) THEN Y(J+1)=Ynew IF ((TJ+H).GT.Br) THEN T(J+1)=B ELSE T(J+1)=TJ+H ENDIF J=J+1 TJ=T(J) ENDIF IF (Err.EQ.0) THEN S=0 ELSE S=0.84*(Tol*H/Err)**0.25 ENDIF IF (S.LT.0.75.AND.H.GT.(2*Hmin)) THEN H=H/2 ENDIF IF (S.GT.1.50.AND.(2*H).LT.Hmax) THEN H=H*2 ENDIF IF (Big.LT.ABS(Y(J)).OR.MaxM.EQ.J) RETURN Mend=J IF (B.GT.T(J)) THEN M=J+1 ELSE M=J ENDIF WRITE(9,1000) T(J),Y(J) GOTO 10 ENDIF RETURN 1000 FORMAT(5X,2F18.9) END SUBROUTINE MESSAGE(Meth,Order,State) INTEGER K,Meth,Order,State REAL Resp WRITE(9,*)' Solution of the differential equation Y` = F +(T,Y)' WRITE(9,*)' ' WRITE(9,*)' with the initial condition Y(A) = Y .' WRITE(9,*)' 0' WRITE(9,*)' ' WRITE(9,*)' A numerical approximation is computed over [A +,B].' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Choose the method of approximation:' WRITE(9,*)' ' WRITE(9,*)' < 1 > Adams-Bashforth-Moulton method' WRITE(9,*)' ' WRITE(9,*)' < 2 > Milne-Simpson`s method' WRITE(9,*)' ' WRITE(9,*)' < 3 > Hamming`s method' WRITE(9,*)' ' WRITE(9,*)' < 4 > Taylor`s method of order N=3' WRITE(9,*)' ' WRITE(9,*)' < 5 > Taylor`s method of order N=4' WRITE(9,*)' ' WRITE(9,*)' < 6 > Runge-Kutta-Fehlberg' WRITE(9,*)' ' WRITE(9,*)' SELECT < 1 - 6 > ? ' Resp=Meth READ(9,*) Resp Meth=INT(Resp) IF ((Meth.LT.1.OR.Meth.GT.6).AND.State.NE.0) THEN Meth=1 ENDIF IF (Meth.EQ.4) THEN Order=3 ENDIF IF (Meth.EQ.5) THEN Order=4 ENDIF RETURN END SUBROUTINE INPUT(Meth) INTEGER Meth REAL Valu CHARACTER Ans*1 DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' You chose ' WRITE(9,*)' ' IF (Meth.EQ.1) THEN WRITE(9,*)' Adams-Bashforth-Moulton method to solve + the D.E.' ELSEIF (Meth.EQ.2) THEN WRITE(9,*)' the Milne-Simpson method to solve + the D.E.' ELSEIF (Meth.EQ.3) THEN WRITE(9,*)' Hamming`s method to solve the D.E.' ELSEIF (Meth.EQ.4) THEN WRITE(9,*)' Taylor`s method of order N=3 to solve + the D.E.' ELSEIF (Meth.EQ.5) THEN WRITE(9,*)' Taylor`s method of order N=4 to solve + the D.E.' ELSEIF (Meth.EQ.6) THEN WRITE(9,*)' the Runge-Kutta-Fehlberg method to solve + the D.E.' ENDIF WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' with the initial condition Y(A) = Y0.' WRITE(9,*)' ' WRITE(9,*)' A numerical approximation is computed over [A,B].' WRITE(9,*)' ' WRITE(9,*)' You must enter the endpoints for the interval,' WRITE(9,*)' ' WRITE(9,*)' the initial condition Y0, and the number of steps + M.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,'(A)') Ans RETURN END SUBROUTINE EPOINTS(A,B,Y0,M,Tol,State) INTEGER I,M,Stat,State REAL A,B,Tol,Y0,Valu C STATUS = (Change,Enter,Done) CHARACTER Ans*1 Stat=1 IF (State.EQ.0) THEN Stat=0 ENDIF 10 CONTINUE DO 20 I=1,15 WRITE(9,*)' ' 20 CONTINUE CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' IF (Stat.EQ.1) THEN WRITE(9,*)' ENTER the left endpoint A = ' READ(9,*) A WRITE(9,*)' ' WRITE(9,*)' ENTER the right endpoint B = ' READ(9,*) B WRITE(9,*)' ' WRITE(9,*)' ENTER initial condition Y(A) = ' READ(9,*) Y0 WRITE(9,*)' ' WRITE(9,*)' ENTER the number of steps M = ' Valu=M READ(9,*) Valu M=INT(Valu) IF (M.LT.1) THEN M=1 ENDIF IF (M.GT.1000) THEN M=1000 ENDIF WRITE(9,*)' ENTER the tolerance Tol = ' Valu=Tol READ(9,*) Valu Tol=Valu ELSE WRITE(9,*)' The left endpoint is A =',A WRITE(9,*)' ' WRITE(9,*)' The right endpoint is B =',B WRITE(9,*)' ' WRITE(9,*)' Initial condition is Y(A) =',Y0 WRITE(9,*)' ' WRITE(9,*)' The number of steps is M = ',M WRITE(9,*)' ' WRITE(9,*)' The tolerance is Tol = ',Tol ENDIF WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Do you want to make a change ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y'.OR.Ans.EQ.'y') THEN Stat=0 DO 30 I=1,15 WRITE(9,*)' ' 30 CONTINUE WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' The current left endpoint is A =',A WRITE(9,*)' ENTER the NEW left endpoint A = ' READ(9,*) A WRITE(9,*)' ' WRITE(9,*)' The current right endpoint is B =',B WRITE(9,*)' ENTER the NEW right endpoint B = ' READ(9,*) B WRITE(9,*)' ' WRITE(9,*)' The current I. C. is Y(A) =',Y0 WRITE(9,*)' ENTER the NEW I. C. Y(A) = ' READ(9,*) Y0 WRITE(9,*)' ' WRITE(9,*)' The current value of M is M = ',M WRITE(9,*)' ENTER the NEW value of M = ' Valu=M READ(9,*) Valu M=INT(Valu) IF (M.LT.1) THEN M=1 ENDIF IF (M.GT.1000) THEN M=1000 ENDIF WRITE(9,*)' The current value of Tol is Tol = ',Tol WRITE(9,*)' ENTER the NEW value of Tol = ' READ(9,*) Valu Tol=Valu ELSE Stat=2 ENDIF IF (Stat.EQ.0.OR.Stat.EQ.1) GOTO 10 RETURN END SUBROUTINE RESULTS(T,Y,M,Mend,Meth) INTEGER I,K,M,Mend,Meth REAL T,Y DIMENSION T(1:1000),Y(1:1000) DO 10 I=1,10 WRITE(9,*)' ' 10 CONTINUE IF (Meth.EQ.1) THEN WRITE(9,*)'Adams-Bashforth-Moulton method was used to solve the & D.E.' ELSEIF (Meth.EQ.2) THEN WRITE(9,*)'The Milne-Simpson method was used to solve the D.E.' ELSEIF (Meth.EQ.3) THEN WRITE(9,*)'Hamming`s method was used to solve the D.E.' ELSEIF (Meth.EQ.4) THEN WRITE(9,*)'Taylor`s method of order N=3 was used to solve the & D.E.' ELSEIF (Meth.EQ.5) THEN WRITE(9,*)'Taylor`s method of order N=4 was used to solve the & D.E.' ELSEIF (Meth.EQ.6) THEN WRITE(9,*)'The Runge-Kutta-Fehlberg method to solve the D.E.' ENDIF WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)'with Y(',T(0),' ) =',Y(0) WRITE(9,*)' ' WRITE(9,*)' K T(K) Y(K)' WRITE(9,*)' ---------------------------------------' DO 20 K=0,Mend WRITE(9,999) K, T(K), Y(K) 20 CONTINUE IF (Mend.LT.M) THEN WRITE(9,*)' ' WRITE(9,*)'The solution points are approaching a pole.' ENDIF 999 FORMAT(I4,4X,4F15.7) PAUSE RETURN END PROGRAM BDDVALUE 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 9.9 (Linear Shooting Method). C Section 9.8, Boundary Value Problems, Page 488 C C Algorithm 9.10 (Finite-Difference Method). C Section 9.9, Finite-Difference Method, Page 496 C PARAMETER(GNmax = 630,MaxM = 101) CHARACTER Ans*1,Ach*1,Mess*80 REAL D,T,X,X1,X2 DIMENSION D(1:4),T(0:MaxM),X(0:MaxM),X1(0:MaxM),X2(0:MaxM) INTEGER FunM,FunT, GNpts, Inum, M, Mend, Meth, Order, Sub REAL A, B, Alpha, Beta INTEGER State,DoMo,Changes,Done,Working,Go,Stop EXTERNAL F1,F2,P,Q,R FunM = 9 Meth = 1 CALL MESSAGE(Meth) CALL INPUT(Meth,FunT,FunM) CALL EPOINTS(A,B,Alpha,Beta,M) IF Meth.EQ.1 THEN CALL LinSho(A,B,Alpha,Beta,M,T,X,FunT) ENDIF IF Meth.EQ.2 THEN CALL FinDif(A,B,Alpha,Beta,M,T,X,FunT) ENDIF CALL RESULTS(FunT,Meth,T,X,M,Mend) WRITE(9,*)' ' READ(9,*) Ans STOP END !End of the Main Program REAL FUNCTION P(T,FunT) INTEGER FunT REAL T IF (FunT.EQ.1) THEN P = 2.0*T/(1.0 + T*T) ELSEIF (FunT.EQ.2) THEN P = -2.0/T ELSEIF (FunT.EQ.3) THEN P = -5.0 ELSEIF (FunT.EQ.4) THEN P = -2.0 ELSEIF (FunT.EQ.5) THEN P = -1.0/T ELSEIF (FunT.EQ.6) THEN P = -2.0/T ELSEIF (FunT.EQ.7) THEN P = -1.0/T ELSEIF (FunT.EQ.8) THEN P = 1.0/T ELSEIF (FunT.EQ.9) THEN P = -1.0/T ENDIF RETURN END REAL FUNCTION Q(T,FunT) INTEGER FunT REAL T IF (FunT.EQ.1) THEN Q = -2.0/(1.0 + T*T) ELSEIF (FunT.EQ.2) THEN Q = 2.0/(T*T) ELSEIF (FunT.EQ.3) THEN Q = -6.0 ELSEIF (FunT.EQ.4) THEN Q = -2.0 ELSEIF (FunT.EQ.5) THEN Q = -(1.0 - 1.0/(4.0*T*T)) ELSEIF (FunT.EQ.6) THEN Q = 2.0/(T*T) ELSEIF (FunT.EQ.7) THEN Q = -(1.0 - 1.0/(4.0*T*T)) ELSEIF (FunT.EQ.8) THEN Q = -1.0/(T*T) ELSEIF (FunT.EQ.9) THEN Q = -16.0/(T*T) ENDIF RETURN END REAL FUNCTION R(T,FunT) INTEGER FunT REAL T IF (FunT.EQ.1) THEN R = 1.0 ELSEIF (FunT.EQ.2) THEN R = 10.0*COS(LOG(T))/(T*T) ELSEIF (FunT.EQ.3) THEN R = T*EXP(-2.0*T) + 3.9*COS(3.0*T) ELSEIF (FunT.EQ.4) THEN R = EXP(-T) + SIN(2.0*T) ELSEIF (FunT.EQ.5) THEN R = 0.0 ELSEIF (FunT.EQ.6) THEN R = SIN(T)/(T*T) ELSEIF (FunT.EQ.7) THEN R = SQRT(T)*COS(T) ELSEIF (FunT.EQ.8) THEN R = 1.0 ELSEIF (FunT.EQ.9) THEN R = 1.0/(T*T) ENDIF RETURN END REAL FUNCTION F1(T,X,Y,FunT) INTEGER FunT REAL T,X,Y IF (FunT.EQ.1) THEN F1 = 2.0*T*Y/(1.0 + T*T) - 2.0*X/(1.0 + T*T) + 1.0 ELSEIF (FunT.EQ.2) THEN F1 = -2.0/T*Y + 2.0/(T*T)*X + 10.0*COS(LOG(T))/(T*T) ELSEIF (FunT.EQ.3) THEN F1 = -5.0*Y - 6.0*X + T*EXP(-2.0*T) + 3.9*COS(3.0*T) ELSEIF (FunT.EQ.4) THEN F1 = -2.0*Y - 2.0*X + EXP(-T) + SIN(2.0*T) ELSEIF (FunT.EQ.5) THEN F1 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X ELSEIF (FunT.EQ.6) THEN F1 = -2.0/T*Y + 2.0/(T*T)*X + SIN(T)/(T*T) ELSEIF (FunT.EQ.7) THEN F1 = -1.0/T*Y-(1.0 - 1.0/(4.0*T*T))*X+SQRT(T)*COS(T) ELSEIF (FunT.EQ.8) THEN F1 = 1.0/T*Y - 1.0/(T*T)*X + 1.0 ELSEIF (FunT.EQ.9) THEN F1 = -1.0/T*Y - 16.0/(T*T)*X + 1.0/(T*T) ENDIF RETURN END REAL FUNCTION F2(T,X,Y,FunT) INTEGER FunT REAL T,X,Y IF (FunT.EQ.1) THEN F2 = 2.0*T*Y/(1.0 + T*T) - 2.0*X/(1.0 + T*T) ELSEIF (FunT.EQ.2) THEN F2 = -2.0/T*Y + 2.0/(T*T)*X ELSEIF (FunT.EQ.3) THEN F2 = -5.0*Y - 6.0*X ELSEIF (FunT.EQ.4) THEN F2 = -2.0*Y - 2.0*X ELSEIF (FunT.EQ.5) THEN F2 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X ELSEIF (FunT.EQ.6) THEN F2 = -2.0/T*Y + 2.0/(T*T)*X ELSEIF (FunT.EQ.7) THEN F2 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X ELSEIF (FunT.EQ.8) THEN F2 = 1.0/T*Y - 1.0/(T*T)*X ELSEIF (FunT.EQ.9) THEN F2 = -1.0/T*Y - 16/(T*T)*X ENDIF RETURN END SUBROUTINE PRTFUN(FunT) INTEGER FunT IF (FunT.EQ.1) THEN WRITE(9,*)'X`` = 2*T*X`/(1 + T*T) - 2*X/(1 + T*T) + 1' ELSEIF (FunT.EQ.2) THEN WRITE(9,*)'X`` = -2/T*X` + 2/(T*T)*X + 10*COS(LOG(T))/(T*T)' ELSEIF (FunT.EQ.3) THEN WRITE(9,*)'X`` = -5*X` - 6*X + T*EXP(-2*T) + 3.9*COS(3*T)' ELSEIF (FunT.EQ.4) THEN WRITE(9,*)'X`` = -2*X` - 2*X + EXP(-T) + SIN(2*T)' ELSEIF (FunT.EQ.5) THEN WRITE(9,*)'X`` = -1/T*X` - (1 - 1/(4*T*T))*X' ELSEIF (FunT.EQ.6) THEN WRITE(9,*)'X`` = -2/T*X` + 2/(T*T)*X + SIN(T)/(T*T)' ELSEIF (FunT.EQ.7) THEN WRITE(9,*)'X`` = -1/T*X` - (1 - 1/(4*T*T))*X + SQRT(T)*COS(T)' ELSEIF (FunT.EQ.8) THEN WRITE(9,*)'X`` = 1/T*X` - 1/(T*T)*X + 1' ELSEIF (FunT.EQ.9) THEN WRITE(9,*)'X`` = -1/T*X` - 16/(T*T)*X + 1/(T*T)' ENDIF RETURN END !Start of the subroutine for the Linear Shooting method. SUBROUTINE RKbdF(F,A,B,Alpha,Beta,M,T,X,FunT) PARAMETER(MaxM = 101) INTEGER J,M,FunT REAL A,B,Alpha,Beta,T,X,Y REAL H,S1,S2 S3,S4,R1,R2,R3,R4,Tj,Xj,Yj DIMENSION T(0:MaxM),X(0:MaxM),Y(0:MaxM) EXTERNAL F H = (B - A)/M T(0) = A X(0) = Alpha Y(0) = Beta DO J=0,(M-1) Tj = T(J) Xj = X(J) Yj = Y(J) S1 = H*Yj R1 = H*F(Tj, Xj, Yj,FunT) S2 = H*(Yj + R1/2.0) R2 = H*F(Tj + H/2.0, Xj+S1/2.0, Yj+R1/2.0,FunT) S3 = H*(Yj + R2/2.0) R3 = H*F(Tj + H/2.0, Xj+S2/2.0, Yj+R2/2.0,FunT) S4 = H*(Yj + R3) R4 = H*F(Tj + H, Xj + S3, Yj + R3,FunT) X(J+1) = Xj + (S1 + 2.0*S2 + 2.0*S3 + S4)/6.0 Y(J+1) = Yj + (R1 + 2.0*R2 + 2.0*R3 + R4)/6.0 T(J+1) = A + H*(J+1) ENDDO END !Start of the main program for the Linear Shooting method. SUBROUTINE LinSho(A,B,Alpha,Beta,M,T,X,FunT) PARAMETER(MaxM = 101) INTEGER J,M,FunT REAL A,B,Alpha,Beta REAL T,X,X1,X2 DIMENSION T(0:MaxM),X(0:MaxM),X1(0:MaxM),X2(0:MaxM) EXTERNAL F1,F2 CALL RKbdF(F1,A,B,Alpha,0.0,M,T,X1,FunT) CALL RKbdF(F2,A,B, 0.0 ,1.0,M,T,X2,FunT) DO J=0,M X(J) = X1(J) + (Beta - X1(M))*X2(J)/X2(M) ENDDO END !The main program for the Linear Shooting method. !Start of the subroutines for the Finite Difference method. SUBROUTINE CoeMat(A,B,Alpha,Beta,N,Va,Vb,Vc,Vd,Vt,FunT) PARAMETER(MaxM = 101) INTEGER J,N,FunT REAL A,B,Alpha,Beta REAL Va, Vb, Vc, Vd, Vt DIMENSION Va(0:MaxM),Vb(0:MaxM),Vc(0:MaxM),Vd(0:MaxM),Vt(0:MaxM) EXTERNAL P,Q,R H = (B - A)/N DO J=1,(N-1) Vt(J) = A + H*J ENDDO DO J=1,(N-1) Vb(J) = -H*H*R(Vt(J),FunT) ENDDO Vb(1) = Vb(1) + (1 + H/2*P(Vt(1),FunT))*Alpha Vb(N-1) = Vb(N-1) + (1 - H/2*P(Vt(N-1),FunT))*Beta DO J=1,(N-1) Vd(J) = 2 + H*H*Q(Vt(J),FunT) ENDDO DO J=1,(N-2) Va(J) = -1 - H/2*P(Vt(J+1),FunT) ENDDO DO J=1,(N-2) Vc(J) = -1 + H/2*P(Vt(J),FunT) ENDDO END SUBROUTINE TriMat(Va,Vb,Vc,Vd,N,X) PARAMETER(MaxM = 101) INTEGER K,N REAL Va,Vb,Vc,Vd,X DIMENSION Va(0:MaxM),Vb(0:MaxM),Vc(0:MaxM),Vd(0:MaxM),X(0:MaxM) REAL TEM DO K=2,(N-1) TEM = Va(K-1)/Vd(K-1) Vd(K) = Vd(K) - TEM*Vc(K-1) Vb(K) = Vb(K) - TEM*Vb(K-1) ENDDO X(N-1) = Vb(N-1)/Vd(N-1) DO K=(N-2),1,-1 X(K) = (Vb(K) - Vc(K)*X(K+1))/Vd(K) ENDDO END !Start of the main program for the Finite Difference method. SUBROUTINE FinDif(A,B,Alpha,Beta,N,T,X,FunT) PARAMETER(MaxM = 101) INTEGER N,FunT REAL A,B,Alpha,Beta,H REAL T, X DIMENSION T(0:MaxM), X(0:MaxM) REAL Va, Vb, Vc, Vd DIMENSION Va(0:MaxM), Vb(0:MaxM), Vc(0:MaxM), Vd(0:MaxM) !Remark: Some of the zero components are NOT used. EXTERNAL P,Q,R CALL CoeMat(A,B,Alpha,Beta,N,Va,Vb,Vc,Vd,T,FunT) CALL TriMat(Va, Vb, Vc, Vd, N, X) T(0) = A T(N) = B X(0) = Alpha X(N) = Beta END !The main program for the Finite Difference method. SUBROUTINE MESSAGE(Meth) INTEGER Meth CHARACTER Mess*80 WRITE(9,*)' BOUNDARY VALUE DIFFERENTIAL EQUATIONS' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Solution of the differential equation' WRITE(9,*)' ' WRITE(9,*)'x`` = p(t)x ` ( t ) + q ( t ) x ( t ) + r ( t )' WRITE(9,*)' ' WRITE(9,*)'with x(a)=Alpha and x(b)=Beta computed over [a,b].' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Choose the method of approximation:' WRITE(9,*)' ' WRITE(9,*)' < 1 > Linear shooting method' WRITE(9,*)' ' WRITE(9,*)' < 2 > Finite difference method' WRITE(9,*)' ' Mess = ' SELECT < 1 - 2 > ? ' Meth = 1 WRITE(9,*) Mess READ(9,*) Meth IF (Meth .LE. 1) .AND. (State .NE. Changes) THEN Meth = 1 ENDIF IF (Meth .GT. 2) .AND. (State .NE. Changes) THEN Meth = 2 ENDIF END SUBROUTINE INPUT(Meth,FunT,FunM) INTEGER FunT,FunM,K CHARACTER Ans*1,Mess*80 WRITE(9,*)' ' IF Meth.EQ.1 THEN WRITE(9,*)'You chose the linear shooting method' ENDIF IF Meth.EQ.2 THEN WRITE(9,*)'You chose the finite difference method' ENDIF WRITE(9,*)'to solve Y` = F(T,Y).' WRITE(9,*)' ' WRITE(9,*)'Choose your D.E.:' WRITE(9,*)' ' DO K=1,FunM WRITE(9,*)'<', K,'>' CALL PRTFUN(K) WRITE(9,*)' ' ENDDO Mess = ' SELECT < 1 - 9 > ? ' WRITE(9,*) Mess READ(9,*) FunT IF (FunT.LT.1) THEN FunT = 1 ENDIF IF (FunT.GT.FunM) THEN FunT = FunM ENDIF WRITE(9,*)' ' IF Meth.EQ.1 THEN WRITE(9,*)'You chose the linear shooting method' ENDIF IF Meth.EQ.2 THEN WRITE(9,*)'You chose the finite difference method' ENDIF WRITE(9,*)'to solve the D.E.' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunT) WRITE(9,*)' ' WRITE(9,*)'with the initial conditions' WRITE(9,*)' ' WRITE(9,*)'X(A) = Alpha and X(B) = Beta.' WRITE(9,*)' ' WRITE(9,*)'A numerical approximation is computed over [A,B].' WRITE(9,*)' ' WRITE(9,*)'You must give the endpoints for the interval,' WRITE(9,*)' ' WRITE(9,*)'the boundary values, and the number of steps.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,*) Ans WRITE(9,*)' ' end SUBROUTINE EPOINTS(A,B,Alpha,Beta,M) INTEGER M REAL A,B,Alpha,Beta CHARACTER Mess*80 WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunT) WRITE(9,*)' ' WRITE(9,*)' ' Mess = 'ENTER the left endpoint A = ' WRITE(9,*) Mess READ(9,*) A WRITE(9,*)' ' Mess = 'ENTER the right endpoint B = ' WRITE(9,*) Mess READ(9,*) B WRITE(9,*)' ' Mess = 'ENTER initial condition X(A) = ' WRITE(9,*) Mess READ(9,*) Alpha WRITE(9,*)' ' Mess = 'ENTER initial condition X(B) = ' WRITE(9,*) Mess READ(9,*) Beta WRITE(9,*)' ' Mess = 'ENTER the number of steps M = ' M = 1 WRITE(9,*) Mess READ(9,*) M WRITE(9,*)' ' IF (M.LT.1) M = 1 IF (M.GT.1000) M = 1000 END SUBROUTINE RESULTS(FunT,Meth,T,X,M,Mend) INTEGER FunT,K,M,Mend REAL T,X DIMENSION T(0:MaxM),X(0:MaxM) WRITE(9,*)' ' IF Meth.EQ.1 THEN WRITE(9,*)' The linear shooting method' ENDIF IF Meth.EQ.2 THEN WRITE(9,*)' The finite difference method' ENDIF WRITE(9,*)' was used to solve the D.E.' WRITE(9,*)' ' CALL PRTFUN(FunT) WRITE(9,*)' ' WRITE(9,*)'with X(', T(0), ' ) =', X(0) WRITE(9,*)'and X(', T(M), ' ) =', X(M) WRITE(9,*)' ' WRITE(9,*)' K', ' T(K) ', ' X(K)' WRITE(9,*)' --------------------------------------------' WRITE(9,*)' ' DO K=0,M WRITE(9,*)' ',K, ' ', T(K), ' ', X(K) ENDDO Mend = M IF (Mend.LT.M) THEN WRITE(9,*)'The solution points are approaching a pole.' WRITE(9,*)' ' ENDIF WRITE(9,*)' Press the key. ' READ(9,*) Ans WRITE(9,*)' ' END