PROGRAM NUMDERIV 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 6.1 (Differentiation Using Limits). C Section 6.1, Approximating the Derivative, Page 326 C PARAMETER(Tol=1E-7,Max=100) INTEGER M REAL DX,H,X DIMENSION D(0:Max),DX(0:Max),E(0:Max) CHARACTER ANS*1 EXTERNAL F 10 CALL INPUTS(X) CALL DLIMIT(F,D,DX,E,X,H,M,Tol) CALL RESULT(D,DX,E,X,H,M) WRITE(9,*)' ' WRITE(9,*) 'WANT TO COMPUTE ANOTHER VALUE OF F`(X) ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X) REAL X F=COS(X) RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' F(X) = COS(X)' RETURN END SUBROUTINE DLIMIT(F,D,DX,E,X,H,M,Tol) PARAMETER(Max=100) INTEGER M,N REAL D,DX,E,H,R,X DIMENSION D(0:Max),DX(0:Max),E(0:Max),R(0:Max) CHARACTER ANS*1 EXTERNAL F Small=1E-9 H = 1 DX(0) = H D(0) = 0.5*(F(X+H)-F(X-H))/H DO N=1,2 H = H/2 DX(N) = H D(N) = 0.5*(F(X+H)-F(X-H))/H E(N) = ABS(D(N)-D(N-1)) R(N) = 2*E(N)/(ABS(D(N))+ABS(D(N-1))+Small) ENDDO N = 1 WHILE ((E(N).GT.E(N+1) .OR. R(N).GT.Tol).AND.(N.LT.Max)) H = H/2 DX(N+2) = H D(N+2) = 0.5*(F(X+H)-F(X-H))/H E(N+2) = ABS(D(N+2)-D(N+1)) R(N+2) = 2*E(N+2)/(ABS(D(N+2))+ABS(D(N+1))+Small) N = N+1 REPEAT M=N RETURN END SUBROUTINE XLIMIT(F,D,DX,E,X,H,M,Tol) C This subroutine uses labeled DO loop(s). PARAMETER(Max=100) INTEGER M,N REAL D,DX,E,H,R,X DIMENSION D(0:Max),DX(0:Max),E(0:Max),R(0:Max) CHARACTER ANS*1 EXTERNAL F Small=1E-9 H = 1 DX(0) = H D(0) = 0.5*(F(X+H)-F(X-H))/H DO 10 N=1,2 H = H/2 DX(N) = H D(N) = 0.5*(F(X+H)-F(X-H))/H E(N) = ABS(D(N)-D(N-1)) R(N) = 2*E(N)/(ABS(D(N))+ABS(D(N-1))+Small) 10 CONTINUE N = 1 20 IF ((E(N).GT.E(N+1) .OR. R(N).GT.Tol).AND.(N.LT.Max)) THEN H = H/2 DX(N+2) = H D(N+2) = 0.5*(F(X+H)-F(X-H))/H E(N+2) = ABS(D(N+2)-D(N+1)) R(N+2) = 2*E(N+2)/(ABS(D(N+2))+ABS(D(N+1))+Small) N = N+1 GOTO 20 ENDIF M=N RETURN END SUBROUTINE INPUTS(X) INTEGER I REAL X DO 10 I=1,16 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE LIMIT OF THE SYMMETRIC DIFFERENCE QUOTIENTS' WRITE(9,*)' ' WRITE(9,*)'IS USED TO FIND AN APPROXIMATION TO F`(X):' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' F`(X) = LIM [F(X+H) - F(X-H)]/2H' WRITE(9,*)' H -> 0 ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'THE GIVEN FUNCTION IS:' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE VALUE X = ' READ(9,*) X RETURN END SUBROUTINE RESULT(D,DX,E,X,H,M) PARAMETER(Max=100) INTEGER I,M REAL B,D,DL,DR,DX,E,H,X DIMENSION D(0:Max),DX(0:Max),E(0:Max) EXTERNAL F WRITE(9,*)' ' WRITE(9,*)' H & ( F(X+H) - F(X-H) )/2H' DO 10 N=0,M-1 WRITE(9,1000) DX(N),D(N) 10 CONTINUE WRITE(9,*)' ' PAUSE WRITE(9,*)' ' WRITE(9,*)'THE LIMIT OF THE SYMMETRIC DIFFERENCE QUOTIENTS WAS + FOUND' WRITE(9,*)' ' WRITE(9,*)' F`(X) = LIM [F(X+H) - F(X-H)]/2H' WRITE(9,*)' H -> 0 ' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' AFTER ',M,' ITERATIONS, H =',H WRITE(9,*)' ' WRITE(9,*)'AN APPROXIMATION FOR THE VALUE OF THE DERIVATIVE IS:' WRITE(9,*)' ' WRITE(9,*)'F`(',X,' ) =',D(M-1) WRITE(9,*)' ' WRITE(9,*)' THE ACCURACY IS +-',E(M-1) DL=(F(X)-F(X-H))/H DR=(F(X+H)-F(X))/H B=1+ABS(DL)+ABS(DR) IF (ABS(DR-DL)>0.001*B) THEN WRITE(9,*)' ' WRITE(9,*)'THE LEFT-HAND AND RIGHT-HAND DERIVATIVES ARE NOT + EQUAL.' WRITE(9,*)' ' WRITE(9,*)'LEFT-HAND [F(X)-F(X-H)]/H =',DL WRITE(9,*)' ' WRITE(9,*)'RIGHT-HAND [F(X+H)-F(X)]/H =',DR WRITE(9,*)' ' WRITE(9,*)'PERHAPS F`(X) DOES NOT EXIST.' ENDIF PAUSE RETURN 1000 FORMAT(F20.10,5X,F20.10) END PROGRAM DIFFEXTR 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 6.2 (Differentiation Using Extrapolation). C Section 6.1, Approximating the Derivative, Page 327 C PARAMETER(Tol=1E-7,Delta=1E-7,MaxN=15) INTEGER J,K,N REAL D,DX,Error,H,RelErr,X DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN) CHARACTER ANS*1 EXTERNAL F 10 CALL INPUTS(X) CALL EXTRAP(F,D,DX,X,H,Error,N,Tol,Delta) CALL RESULT(F,D,DX,X,H,Error,N) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE F`(X) AT ANOTHER POINT ? ' 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) RETURN END SUBROUTINE PRINTFUN WRITE(9,*)' F(X) = SIN(X)' RETURN END SUBROUTINE EXTRAP(F,D,DX,X,H,Error,N,Tol,Delta) PARAMETER(MaxN=15) INTEGER J,K,N REAL D,DX,Error,H,RelErr,Small,X DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN) EXTERNAL F Small=1E-7 H=1 DX(0) = H J=1 Error=1 RelErr=1 D(0,0)=0.5*(F(X+H)-F(X-H))/H WHILE ( RelErr.GT.Tol .AND. Error.GT.Delta .AND. J.LT.16) H=H/2 DX(J) = H D(J,0)=0.5*(F(X+H)-F(X-H))/H DO K=1,J D(J,K)=D(J,K-1)+(D(J,K-1)-D(J-1,K-1))/(4**K-1) ENDDO Error=ABS(D(J,J)-D(J-1,J-1)) RelErr=2*Error/(ABS(D(J,J))+ABS(D(J-1,J-1))+Small) N=J J=J+1 REPEAT RETURN END SUBROUTINE XEXTRAP(F,D,DX,X,H,Error,N,Tol,Delta) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL D,DX,Error,H,RelErr,Small,X DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN) EXTERNAL F Small=1E-7 H=1 DX(0) = H J=1 Error=1 RelErr=1 D(0,0)=0.5*(F(X+H)-F(X-H))/H 10 IF ( RelErr.GT.Tol .AND. Error.GT.Delta .AND. J.LT.16) THEN H=H/2 DX(J) = H D(J,0)=0.5*(F(X+H)-F(X-H))/H DO 20 K=1,J D(J,K)=D(J,K-1)+(D(J,K-1)-D(J-1,K-1))/(4**K-1) 20 CONTINUE Error=ABS(D(J,J)-D(J-1,J-1)) RelErr=2*Error/(ABS(D(J,J))+ABS(D(J-1,J-1))+Small) N=J J=J+1 GOTO 10 ENDIF RETURN END SUBROUTINE INPUTS(X) INTEGER I REAL X DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'RICHARDSON`S EXTRAPOLATION APPLIED TO THE SEQUENCE OF + DIFFERENCE QUOTIENTS' WRITE(9,*)' ' WRITE(9,*)' -K ' WRITE(9,*)' { [F(X+H ) - F(X-H )]/2H : WHERE H = 2 }' WRITE(9,*)' K K K K' WRITE(9,*)' ' WRITE(9,*)'IS USED TO FIND AN APPROXIMATION TO THE DERIVATIVE OF + THE FUNCTION:' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE VALUE X = ' READ(9,*) X RETURN END SUBROUTINE RESULT(F,D,DX,X,H,Error,N) PARAMETER(MaxN=15) INTEGER I,N REAL B,D,DL,DR,DX,Error,H,H0,X DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN) EXTERNAL F WRITE(9,*)' ' WRITE(9,*)' H D(N,N)' DO 10 I=0,N WRITE(9,1000) DX(I),D(I,I) 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'RICHARDSON`S EXTRAPOLATION WAS USED ON THE SEQUENCE OF + DIFFERENCE QUOTIENTS' WRITE(9,*)' -K' WRITE(9,*)' { [F(X+H ) - F(X-H )]/2H : WHERE H = 2 +}' WRITE(9,*)' K K K K' WRITE(9,*)'TO FIND AN APPROXIMATION TO THE DERIVATIVE OF THE + FUNCTION:' WRITE(9,*)' ' CALL PRINTFUN WRITE(9,*)' ' WRITE(9,*)' AFTER ',N,' ITERATIONS, H =',H WRITE(9,*)' ' WRITE(9,*)'AN APPROXIMATION FOR THE VALUE OF THE DERIVATIVE IS' WRITE(9,*)' ' WRITE(9,*)'F`(',X,' ) =',D(N,N) WRITE(9,*)' ' WRITE(9,*)' THE ACCURACY IS +-',Error H0=0.002 IF (H.LT.H0) THEN H0=H ENDIF DL=(F(X)-F(X-H0))/H0 DR=(F(X+H0)-F(X))/H0 B=1+ABS(DL)+ABS(DR) IF (ABS(DR-DL).GT.0.02*B) THEN WRITE(9,*)' ' WRITE(9,*)'THE LEFT-HAND AND RIGHT-HAND DERIVATIVES ARE NOT + EQUAL.' WRITE(9,*)' ' WRITE(9,*)'LEFT-HAND [F(X)-F(X-H0)]/H0 =',DL WRITE(9,*)' ' WRITE(9,*)'RIGHT-HAND [F(X+H0)-F(X)]/H0 =',DR WRITE(9,*)' ' WRITE(9,*)'PERHAPS F`(X) DOES NOT EXIST.' ENDIF PAUSE RETURN 1000 FORMAT(F20.10,5X,F20.10) END PROGRAM NEWPOLYD 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 6.3 (Differentiation Based on N+1 Nodes). C Section 6.2, Numerical Differentiation Formulas, Page 342 C PARAMETER(MaxN=15) INTEGER N REAL A,Df,X,Y CHARACTER ANS*1 DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) 10 CALL GETDATA(X,Y,N) CALL DIFPOLY(A,X,Y,N,Df) CALL RESULTS(A,X,Y,N,Df) WRITE(9,*)' ' WRITE(9,*)'WANT TO USE A DIFFERENT SET OF DATA POINTS ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 END SUBROUTINE DIFPOLY(A,X,Y,N,Df) PARAMETER(MaxN=15) INTEGER J,K,N REAL A,Df,Prod,T,X,Y DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) DO K=0,N A(K)=Y(K) ENDDO DO J=1,N DO K=N,J,-1 A(K)=(A(K)-A(K-1))/(X(K)-X(K-J)) ENDDO ENDDO T=X(0) Df=A(1) Prod=1 DO K=2,N Prod=Prod*(T-X(K-1)) Df=Df+Prod*A(K) ENDDO RETURN END SUBROUTINE XDIFPOLY(A,X,Y,N,Df) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL A,Df,Prod,T,X,Y DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) DO 10 K=0,N A(K)=Y(K) 10 CONTINUE DO 30 J=1,N DO 20 K=N,J,-1 A(K)=(A(K)-A(K-1))/(X(K)-X(K-J)) 20 CONTINUE 30 CONTINUE T=X(0) Df=A(1) Prod=1 DO 40 K=2,N Prod=Prod*(T-X(K-1)) Df=Df+Prod*A(K) 40 CONTINUE RETURN END SUBROUTINE PRINTPOL(A,X,Y,N) PARAMETER(MaxN=15) INTEGER K,N REAL A,X,Y DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) WRITE(9,*)' ' IF (N.EQ.1) THEN WRITE(9,*)'P(T) = A + A [T-X ]' WRITE(9,*)' 0 1 0 ' ELSEIF (N.EQ.2) THEN WRITE(9,*)'P(T) = A + A [T-X ] + A [T-X ][T-X ]' WRITE(9,*)' 0 1 0 2 0 1 ' ELSEIF (N.EQ.3) THEN WRITE(9,*)'P(T) = A + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X +][T-X ]' WRITE(9,*)' 0 1 0 2 0 1 3 0 1 + 2 ' ELSE WRITE(9,*)'P(T) = A + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X +][T-X ]' WRITE(9,*)' 0 1 0 2 0 1 3 0 1 + 2' WRITE(9,*)' +...+ A [T-X ][T-X ] ... [T-X ]' WRITE(9,*)' ',N,' 0 1',(N-1) ENDIF WRITE(9,*)' ' WRITE(9,*)' THE COEFFICIENTS: THE ABSCISSAS: THE OR +DINATES:' WRITE(9,*)'K A(K) X(K) Y( +K)' WRITE(9,*)' ' DO 10 K=0,N WRITE(9,999) K,A(K),X(K),Y(K) 10 CONTINUE 999 FORMAT(1X,I2,3F20.7) RETURN END SUBROUTINE GETDATA(X,Y,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL X,Y DIMENSION X(0:MaxN),Y(0:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' THE DERIVATIVE OF THE NEWTON POLYNOMIAL' WRITE(9,*)' ' WRITE(9,*)'P(T) = A + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X ][ +T-X ]' WRITE(9,*)' 0 1 0 2 0 1 3 0 1 + 2' WRITE(9,*)' +...+ A [T-X ][T-X ]...[T-X ]' WRITE(9,*)' N 0 1 N-1' WRITE(9,*)' ' WRITE(9,*)'IS USED TO APPROXIMATE THE DERIVATIVE AT THE POINT (X , +Y ).' WRITE(9,*)' 0 + 0' WRITE(9,*)' ' WRITE(9,*)'P(T) IS BASED ON THE N+1 NODES (X ,Y ),(X ,Y ),...,(X , +Y ).' WRITE(9,*)' 0 0 1 1 N + N ' WRITE(9,*)' ' WRITE(9,*)'THE VALUES {X } DO NOT NEED TO BE INCREASING.' WRITE(9,*)' K' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',N+1,' POINTS' WRITE(9,*)' ' DO 20 K=0,N WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) WRITE(9,*)'Y(',K,') = ' READ(9,*) Y(K) WRITE(9,*)' ' 20 CONTINUE RETURN END SUBROUTINE RESULTS(A,X,Y,N,Df) PARAMETER(MaxN=15) INTEGER I,N REAL A,Df,X,Y DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' THE NEWTON POLYNOMIAL OF DEGREE ',N,' IS:' WRITE(9,*)' ' CALL PRINTPOL(A,X,Y,N) WRITE(9,*)' ' WRITE(9,*)'THE APPROXIMATION FOR THE DERIVATIVE AT X =',X(0) WRITE(9,*)' 0' WRITE(9,*)'USING THE NEWTON POLYNOMIAL IS:' WRITE(9,*)' ' WRITE(9,*)'P`(',X(0),' ) = ',Df WRITE(9,*)' ' RETURN END