PROGRAM TSERIES 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 4.1 (Evaluation of a Taylor Series). C Section 4.1, Taylor Series and Calculation of Functions, Page 203 C PARAMETER (MAX=501,Tol=0.00000001) INTEGER K,METH,N REAL Close,Sum,X,X0 REAL D(0:MAX) CHARACTER ANS*1 CALL INX0XN(X0,N) CALL GETDERIV(N,METH,D) 10 CALL INPUTX(X) CALL TAYLOR(X0,D,X,Tol,Sum,Close,K,N) CALL RESULTSS(X0,X,Tol,Sum,Close,K,N) WRITE(9,*)' ' WRITE(9,*)'WANT TO TRY ANOTHER X ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END SUBROUTINE TAYLOR(X0,D,X,Tol,Sum,Close,K,N) PARAMETER (MAX=501) INTEGER K,N REAL Close,Prod,Sum,Term,Tol,X,X0 REAL D(0:MAX) Close=1 K=0 Sum=D(0) Prod=1 IF (X.EQ.X0) Close=0 WHILE ((Close.GE.Tol).AND.(K.LT.N)) K=K+1 Prod=Prod*(X-X0)/K WHILE ((D(K).EQ.0).AND.(K.LT.N)) K=K+1 Prod=Prod*(X-X0)/K REPEAT Term=D(K)*Prod IF (Term.NE.0) Close=ABS(Term) Sum=Sum+Term REPEAT RETURN END SUBROUTINE XTAYLOR(X0,D,X,Tol,Sum,Close,K,N) C This subroutine uses simulated WHILE loop(s). PARAMETER (MAX=501) INTEGER K,N REAL Close,Prod,Sum,Term,Tol,X,X0 REAL D(0:MAX) Close=1 K=0 Sum=D(0) Prod=1 IF (X.EQ.X0) Close=0 10 IF ((Close.GE.Tol).AND.(K.LT.N)) THEN K=K+1 Prod=Prod*(X-X0)/K 20 IF ((D(K).EQ.0).AND.(K.LT.N)) THEN K=K+1 Prod=Prod*(X-X0)/K GOTO 20 ENDIF Term=D(K)*Prod IF (Term.NE.0) Close=ABS(Term) Sum=Sum+Term GOTO 10 ENDIF RETURN END SUBROUTINE GETDERIV(N,METH,D) PARAMETER (MAX=501) INTEGER I,K,METH,N REAL D(0:MAX) DO 5 I=1,14 WRITE(9,*)' ' 5 CONTINUE WRITE(9,*)'CHOOSE ONE OF THE FOLLOWING:' WRITE(9,*)' ' WRITE(9,*)'<1> F(X) = EXP(X) EXPANDED ABOUT X0 = 0' WRITE(9,*)' ' WRITE(9,*)'<2> F(X) = COS(X) EXPANDED ABOUT X0 = 0' WRITE(9,*)' ' WRITE(9,*)'<3> F(X) = SIN(X) EXPANDED ABOUT X0 = 0' WRITE(9,*)' ' WRITE(9,*)'<4> ENTER THE DERIVATIVES OF ' WRITE(9,*)' ' WRITE(9,*)' F(X) EXPANDED ABOUT X0 = ',X0 WRITE(9,*)' ' WRITE(9,*)'SELECT <1-4> ? ' READ(9,*) METH WRITE(9,*)' ' GOTO (10,20,30,40), METH 10 CONTINUE C { EXP(X) } X0=0 DO 15 K=0,N D(K)=1 15 CONTINUE GOTO 50 20 CONTINUE C { COS(X) } X0=0 DO 25 K=0,INT(N/4) D(4*K)=1 D(4*K+1)=0 D(4*K+2)=-1 D(4*K+3)=0 25 CONTINUE GOTO 50 30 CONTINUE C { SIN(X) } X0=0 DO 35 K=0,INT(N/4) D(4*K)=0 D(4*K+1)=1 D(4*K+2)=0 D(4*K+3)=-1 35 CONTINUE GOTO 50 40 CONTINUE WRITE(9,*)'ENTER THE ',N+1,' DERIVATIVES OF F(X)' WRITE(9,*)'EVALUATED AT THE CENTER VALUE X0 = ',X0 WRITE(9,*)' ' WRITE(9,*)'D(0) , D(1) , D(2) ,..., D(N)' WRITE(9,*)' ' DO 45 K=0,N WRITE(9,*)'D(',K,') = ' READ(9,*) D(K) 45 CONTINUE 50 CONTINUE RETURN END SUBROUTINE INX0XN(X0,N) PARAMETER (MAX=501) INTEGER I,N REAL X0 DO 10 I=1,14 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'COMPUTER EVALUATION OF THE TAYLOR`S SERIES' WRITE(9,*)' ' WRITE(9,*)'FOR THE FUNCTION F(X) CENTERED AT X ' WRITE(9,*)' 0 ' WRITE(9,*)' ' WRITE(9,*)'THE N-THE PARTIAL SUM IS COMPUTED: ' WRITE(9,*)' ' WRITE(9,*)' 2 K +N ' WRITE(9,*)'P(X) = A +A (X-X )+A (X-X ) +...+A (X-X ) +...+A (X-X ) + ' WRITE(9,*)' 0 1 0 2 0 K 0 N 0 + ' WRITE(9,*)' ' WRITE(9,*)' (K) ' WRITE(9,*)'THE K-COEFFICIENT IS A = F (X )/K! = D(K)/K! ' WRITE(9,*)' K 0 ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE CENTER X OF THE SERIES' WRITE(9,*)' 0 ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE CENTER X0 = ' READ(9,*) X0 WRITE(9,*)' ' WRITE(9,*)'ENTER THE DEGREE N = ' READ(9,*) N IF (N.GT.MAX) THEN N=MAX ENDIF WRITE(9,*)' ' RETURN END SUBROUTINE INPUTX(X) REAL X WRITE(9,*)' ' WRITE(9,*)'ENTER THE VALUE OF THE INDEPENDENT VARIABLE' WRITE(9,*)' ' WRITE(9,*)'X = ' READ(9,*) X RETURN END SUBROUTINE RESULTSS(X0,X,Tol,Sum,Close,K,N) INTEGER I,K,N REAL Close,Sum,Tol,X,X0 DO 5 I=1,12 WRITE(9,*)' ' 5 CONTINUE WRITE(9,*)'YOU CHOSE TO APPROXIMATE ' GOTO (10,20,30,40), METH 10 WRITE(9,*)' EXP(',X,' ).' GOTO 50 20 WRITE(9,*)' COS(',X,' ).' GOTO 50 30 WRITE(9,*)' SIN(',X,' ).' GOTO 50 40 WRITE(9,*)'YOUR OWN TAYLOR SERIES.' WRITE(9,*)'EXPANDED ABOUT THE VALUE X0 = ',X0 50 CONTINUE WRITE(9,*)' ' IF ((Close.LE.Tol).AND.(K.LE.N)) THEN WRITE(9,*)'THE VALUE OF THE TAYLOR POLYNOMIAL IS:' WRITE(9,*)' ' WRITE(9,*)' P (',X,' ) =',Sum WRITE(9,*)' ', K WRITE(9,*)' ' WRITE(9,*)'CONSECUTIVE PARTIAL SUMS ARE CLOSER THAN',Close WRITE(9,*)' ' ELSE WRITE(9,*)'THE CURRENT PARTIAL SUM IS:' WRITE(9,*)' ' WRITE(9,*)' P (',X,' ) =',Sum WRITE(9,*)' ',K WRITE(9,*)' ' WRITE(9,*)'CONSECUTIVE PARTIAL SUMS DIFFER BY ',Close WRITE(9,*)' ' WRITE(9,*)'CONVERGENCE HAS NOT BEEN ACHIEVED.' WRITE(9,*)' ' ENDIF RETURN END PROGRAM POLCAL 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 4.2 (Polynomial Calculus). C Section 4.2, Introduction to Interpolation, Page 212 C PARAMETER(MaxN=15) INTEGER Meth,N REAL C,DP,IP,P CHARACTER ANS*1 DIMENSION C(0:MaxN) EXTERNAL DP,IP,P 10 CALL GETPOLY(C,N) 20 CALL CHOICE(C,N,Meth) CALL EVALUATE(C,N,Meth) WRITE(9,*)' ' WRITE(9,*)'WANT ANOTHER CHOICE ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT ANOTHER POLYNOMIAL ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 RETURN END REAL FUNCTION P(C,N,X) PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Poly,X DIMENSION C(0:MaxN) Poly=C(N) DO K=(N-1),0,-1 Poly=C(K)+Poly*X ENDDO P=Poly RETURN END REAL FUNCTION DP(C,N,X) PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Deriv,X DIMENSION C(0:MaxN) Deriv=N*C(N) DO K=(N-1),1,-1 Deriv=K*C(K)+Deriv*X ENDDO DP=Deriv RETURN END REAL FUNCTION IP(C,N,X) PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Integ,X DIMENSION C(0:MaxN) Integ=C(N)/(N+1) DO K=N,1,-1 Integ=C(K-1)/K+Integ*X ENDDO Integ=Integ*X IP=Integ RETURN END REAL FUNCTION XP(C,N,X) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Poly,X DIMENSION C(0:MaxN) Poly=C(N) DO 10 K=(N-1),0,-1 Poly=C(K)+Poly*X 10 CONTINUE P=Poly RETURN END REAL FUNCTION XDP(C,N,X) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Deriv,X DIMENSION C(0:MaxN) Deriv=N*C(N) DO 10 K=(N-1),1,-1 Deriv=K*C(K)+Deriv*X 10 CONTINUE DP=Deriv RETURN END REAL FUNCTION XIP(C,N,X) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER I,K,N REAL C,Integ,X DIMENSION C(0:MaxN) Integ=C(N)/(N+1) DO 10 K=N,1,-1 Integ=C(K-1)/K+Integ*X 10 CONTINUE Integ=Integ*X IP=Integ RETURN END SUBROUTINE PRINTPOL(C,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL C DIMENSION C(0:MaxN) IF (N.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'P(X) = C X + C' WRITE(9,*)' 1 0' WRITE(9,*)' ' ELSEIF (N.EQ.2) THEN WRITE(9,*)' 2' WRITE(9,*)'P(X) = C X + C X + C' WRITE(9,*)' 2 1 0' WRITE(9,*)' ' ELSEIF (N.EQ.3) THEN WRITE(9,*)' 3 2' WRITE(9,*)'P(X) = C X + C X + C X + C ' WRITE(9,*)' 3 2 1 0' WRITE(9,*)' ' ELSE WRITE(9,*)' ',N,' 2' WRITE(9,*)'P(X) = C X +...+ C X + C X + C ' WRITE(9,*)' ',N,' 2 1 0' WRITE(9,*)' ' ENDIF DO 10 K=N,0,-2 IF (K.GE.1) THEN WRITE(9,*)'C(',K,') =',C(K),' C(',(K-1),') =',C(K-1) ELSE WRITE(9,*)'C(',K,') =',C(K) ENDIF 10 CONTINUE RETURN END SUBROUTINE GETPOLY(C,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL C DIMENSION C(0:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'POLYNOMIAL CALCULUS FOR EVALUATING, DIFFERENTIATING' WRITE(9,*)' ' WRITE(9,*)'AND INTEGRATING THE POLYNOMIAL P(X) OF DEGREE N. ' WRITE(9,*)' ' WRITE(9,*)' N N-1 2 ' WRITE(9,*)'P(X) = C X + C X +...+ C X + C X + C ' WRITE(9,*)' N N-1 2 1 0' WRITE(9,*)' ' WRITE(9,*)'ENTER THE DEGREE N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',N+1,' COEFFICIENTS' WRITE(9,*)' ' WRITE(9,*)' C , ... , C , C , C ' WRITE(9,*)' ',N,' 2 1 0' WRITE(9,*)' ' DO 20 K=N,0,-1 WRITE(9,*)'C(',K,') = ' READ(9,*) C(K) 20 CONTINUE RETURN END SUBROUTINE CHOICE(C,N,Meth) PARAMETER(MaxN=15) INTEGER I,Meth,N REAL C DIMENSION C(0:MaxN) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE POLYNOMIAL OF DEGREE ',N,' IS:' WRITE(9,*)' ' CALL PRINTPOL(C,N) WRITE(9,*)' ' WRITE(9,*)'CHOOSE ONE OF THE FOLLOWING:' WRITE(9,*)' ' WRITE(9,*)'1 - EVALUATE P(X)' WRITE(9,*)' ' WRITE(9,*)'2 - EVALUATE P`(X)' WRITE(9,*)' ' WRITE(9,*)'3 - INTEGRATE P(X) OVER [A,B]' WRITE(9,*)' ' WRITE(9,*)'SELECT <1-3> ' READ(9,*) Meth WRITE(9,*)' ' RETURN END SUBROUTINE EVALUATE(C,N,Meth) PARAMETER(MaxN=15) INTEGER I,Meth,N REAL A,B,C,Deriv,DP,IntAB,IP,P,Poly,X DIMENSION C(0:MaxN) EXTERNAL DP,IP,P DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE IF (Meth.EQ.1) THEN WRITE(9,*)'YOU CHOSE TO EVALUATE THE POLYNOMIAL P(X):' ELSEIF (Meth.EQ.2) THEN WRITE(9,*)'YOU CHOSE TO DIFFERENTIATE THE POLYNOMIAL P(X):' ELSEIF (Meth.EQ.3) THEN WRITE(9,*)'YOU CHOSE TO INTEGRATE THE POLYNOMIAL P(X):' ENDIF WRITE(9,*)' ' CALL PRINTPOL(C,N) WRITE(9,*)' ' IF (Meth.EQ.1) THEN WRITE(9,*)'EVALUATE P(X)' WRITE(9,*)' ' WRITE(9,*)' ENTER X = ' READ(9,*) X WRITE(9,*)' ' Poly=P(C,N,X) WRITE(9,*)'P(',X,' ) = ',Poly ELSEIF (Meth.EQ.2) THEN WRITE(9,*)'EVALUATE P`(X)' WRITE(9,*)' ' WRITE(9,*)' ENTER X = ' READ(9,*) X WRITE(9,*)' ' Deriv=DP(C,N,X) WRITE(9,*)'P`(',X,' ) = ',Deriv ELSEIF (Meth.EQ.3) THEN WRITE(9,*)'INTEGRATE P(X) OVER [A,B]' WRITE(9,*)' ' WRITE(9,*)'ENTER THE LEFT ENDPOINT A = ' READ(9,*) A WRITE(9,*)'ENTER THE RIGHT ENDPOINT B = ' READ(9,*) B IntAB=IP(C,N,B)-IP(C,N,A) WRITE(9,*)' ' WRITE(9,*)'B' WRITE(9,*)'/' WRITE(9,*)'| P(X) DX = ',IntAB WRITE(9,*)'/' WRITE(9,*)'A' ENDIF WRITE(9,*)' ' RETURN END PROGRAM LAGPOL 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 4.3 (Lagrange Approximation). C Section 4.3, Lagrange Approximation, Page 224 C PARAMETER(MaxN=15) INTEGER N REAL P,T,X,Y CHARACTER ANS*1 DIMENSION X(0:MaxN),Y(0:MaxN) EXTERNAL P 10 CALL GETDATAS(X,Y,N) 20 CALL EVALUATE(X,Y,N) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE P(T) AGAIN ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO TRY A DIFFERENT SET OF POINTS ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 RETURN END REAL FUNCTION P(X,Y,N,T) PARAMETER(MaxN=15) INTEGER J,K,N REAL Sum,T,Term,X,Y DIMENSION X(0:MaxN),Y(0:MaxN) Sum=0 DO K=0,N Term=Y(K) DO J=0,N IF (J.NE.K) Term=Term*(T-X(J))/(X(K)-X(J)) ENDDO Sum=Sum+Term ENDDO P=Sum RETURN END REAL FUNCTION XP(X,Y,N,T) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL Sum,T,Term,X,Y DIMENSION X(0:MaxN),Y(0:MaxN) Sum=0 DO 20 K=0,N Term=Y(K) DO 10 J=0,N IF (J.NE.K) Term=Term*(T-X(J))/(X(K)-X(J)) 10 CONTINUE Sum=Sum+Term 20 CONTINUE P=Sum RETURN END SUBROUTINE GETDATAS(X,Y,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL X,Y DIMENSION X(0:MaxN),Y(0:MaxN) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'THE LAGRANGE POLYNOMIAL IS USED TO INTERPOLATE' WRITE(9,*)' ' WRITE(9,*)'BASED ON THE N+1 POINTS (X ,Y ),...,(X ,Y )' WRITE(9,*)' 0 0 N N' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE N =',N+1,' POINTS' WRITE(9,*)' ' WRITE(9,*)'(X ,Y ),(X ,Y ),...,(X ,Y )' WRITE(9,*)' 0 0 1 1 N N' 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 EVALUATE(X,Y,N) PARAMETER(MaxN=15) INTEGER K,N REAL T,X,Y DIMENSION X(0:MaxN),Y(0:MaxN) EXTERNAL P DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'THE LAGRANGE POLYNOMIAL INTERPOLATION BASED ON THE N+1 + POINTS' WRITE(9,*)' ' DO 20 K=0,N WRITE(9,*)'X(',K,') =',X(K),' Y(',K,') =',Y(K) 20 CONTINUE WRITE(9,*)' ' WRITE(9,*)'NOW EVALUATE P(T)' WRITE(9,*)' ' WRITE(9,*)'ENTER A VALUE T = ' READ(9,*) T WRITE(9,*)' ' WRITE(9,*)'THE VALUE OF THE LAGRANGE POLYNOMIAL IS ' WRITE(9,*)' ' WRITE(9,*)'P(',T,' ) = ',P(X,Y,N,T) WRITE(9,*)' ' RETURN END PROGRAM NESTEDM 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 4.4 (Nested Multiplication with Multiple Centers). C Section 4.4, Newton Polynomials, Page 233 C PARAMETER(MaxN=15) INTEGER N REAL A,P,T,X CHARACTER ANS*1 DIMENSION A(0:MaxN),X(0:MaxN) EXTERNAL P 10 CALL GETDATAS(A,X,N) 20 CALL EVALUATE(A,X,N,T) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE P(T) AGAIN ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE ANOTHER POLYNOMIAL ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 RETURN END REAL FUNCTION P(A,X,N,T) PARAMETER(MaxN=15) INTEGER K,N REAL A,Sum,T,X DIMENSION A(0:MaxN),X(0:MaxN) Sum=A(N) DO K=(N-1),0,-1 Sum=Sum*(T-X(K))+A(K) ENDDO P= Sum RETURN END REAL FUNCTION XP(A,X,N,T) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER K,N REAL A,Sum,T,X DIMENSION A(0:MaxN),X(0:MaxN) Sum=A(N) DO 10 K=(N-1),0,-1 Sum=Sum*(T-X(K))+A(K) 10 CONTINUE P= Sum RETURN END SUBROUTINE PRINTPOL(A,X,N) PARAMETER(MaxN=15) INTEGER K,N REAL A,X DIMENSION A(0:MaxN),X(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) WRITE(9,*)' ' WRITE(9,*)'N =',N ENDIF WRITE(9,*)' ' WRITE(9,*)'THE COEFFICIENTS: THE CENTERS:' WRITE(9,*)' ' DO 10 K=0,(N-1) WRITE(9,*)'A(',K,') =',A(K),' X(',K,') =', +X(K) 10 CONTINUE WRITE(9,*)'A(',N,') =',A(N) RETURN END SUBROUTINE GETDATAS(A,X,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL A,X DIMENSION A(0:MaxN),X(0:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'NESTED MULTIPLICATION IS USED TO EVALUATE THE NEWTON PO +LYNOMIAL:' 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,*)'THE CENTERS ARE X , X ,..., X AND THE' WRITE(9,*)' 0 1 N-1 ' WRITE(9,*)' ' WRITE(9,*)'COEFFICIENTS ARE A , A ,..., A , A ' WRITE(9,*)' 0 1 N-1 N ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE DEGREE N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',N,' CENTERS' WRITE(9,*)' ' DO 20 K=0,(N-1) WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) 20 CONTINUE WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',N+1,' COEFFICIENTS' WRITE(9,*)' ' DO 30 K=0,N WRITE(9,*)'A(',K,') = ' READ(9,*) A(K) 30 CONTINUE WRITE(9,*)' ' RETURN END SUBROUTINE EVALUATE(A,X,N,T) PARAMETER(MaxN=15) INTEGER I,N REAL A,P,Poly,T,X DIMENSION A(0:MaxN),X(0:MaxN) EXTERNAL P DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE CALL PRINTPOL(A,X,N) WRITE(9,*)' ' WRITE(9,*)'NOW EVALUATE P(T)' WRITE(9,*)' ' WRITE(9,*)'ENTER A VALUE T = ' READ(9,*) T WRITE(9,*)' ' WRITE(9,*)'THE VALUE OF THE NEWTON POLYNOMIAL IS ' WRITE(9,*)' ' Poly=P(A,X,N,T) WRITE(9,*)'P(',T,' ) = ',Poly RETURN END PROGRAM NEPOLY 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 4.5 (Newton Interpolation Polynomial). C Section 4.4, Newton Polynomials, Page 234 C PARAMETER(MaxN=15) INTEGER N REAL A,P,T,X,Y CHARACTER ANS*1 DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) EXTERNAL P 10 CALL GETDATAS(X,Y,N) CALL DIVDIFFS(X,Y,A,N) 20 CALL EVALUATE(X,Y,A,N,T) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE P(T) AGAIN ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE ANOTHER POLYNOMIAL ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END SUBROUTINE DIVDIFFS(X,Y,A,N) PARAMETER(MaxN=15) INTEGER J,K,N REAL A,D,X,Y DIMENSION A(0:MaxN),D(0:MaxN,0:MaxN),X(0:MaxN),Y(0:MaxN) DO K=0,N D(K,0)=Y(K) ENDDO DO J=1,N DO K=J,N D(K,J)=(D(K,J-1)-D(K-1,J-1))/(X(K)-X(K-J)) ENDDO ENDDO DO K=0,N A(K)=D(K,K) ENDDO RETURN END REAL FUNCTION P(A,X,N,T) PARAMETER(MaxN=15) INTEGER K,N REAL A,Sum,T,X DIMENSION A(0:MaxN),X(0:MaxN) Sum=A(N) DO K=(N-1),0,-1 Sum=Sum*(T-X(K))+A(K) ENDDO P=Sum RETURN END SUBROUTINE XDIVDIFFS(X,Y,A,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL A,D,X,Y DIMENSION A(0:MaxN),D(0:MaxN,0:MaxN),X(0:MaxN),Y(0:MaxN) DO 10 K=0,N D(K,0)=Y(K) 10 CONTINUE DO 30 J=1,N DO 20 K=J,N D(K,J)=(D(K,J-1)-D(K-1,J-1))/(X(K)-X(K-J)) 20 CONTINUE 30 CONTINUE DO 40 K=0,N A(K)=D(K,K) 40 CONTINUE RETURN END REAL FUNCTION XP(A,X,N,T) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER K,N REAL A,Sum,T,X DIMENSION A(0:MaxN),X(0:MaxN) Sum=A(N) DO 10 K=(N-1),0,-1 Sum=Sum*(T-X(K))+A(K) 10 CONTINUE P=Sum RETURN END SUBROUTINE GETDATAS(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 NEWTON POLYNOMIAL IS CONSTRUCTED:' 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,*)'BASED ON THE N+1 POINTS (X ,Y ),(X ,Y ),...,(X ,Y )' WRITE(9,*)' 0 0 1 1 N N' WRITE(9,*)' ' WRITE(9,*)'THE CENTERS ARE THE ABSCISSAS X , X ,..., X .' WRITE(9,*)' 0 1 N-1 ' WRITE(9,*)' ' WRITE(9,*)'THE COEFFICIENTS ARE A , A ,..., A , A , AND THEY' WRITE(9,*)' 0 1 N-1 N ' WRITE(9,*)' ' WRITE(9,*)'ARE FOUND BY CONSTRUCTING A DIVIDED DIFFERENCE TABLE.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE DEGREE 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 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 ORDINATES:' 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 EVALUATE(X,Y,A,N,T) PARAMETER(MaxN=15) INTEGER N REAL A,P,T,X DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN) EXTERNAL P DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE CALL PRINTPOL(A,X,Y,N) WRITE(9,*)' ' WRITE(9,*)'NOW EVALUATE P(T)' WRITE(9,*)' ' WRITE(9,*)'ENTER A VALUE T = ' READ(9,*) T WRITE(9,*)' ' WRITE(9,*)'THE VALUE OF THE NEWTON POLYNOMIAL IS ' WRITE(9,*)' ' Poly=P(A,X,N,T) WRITE(9,*)'P(',T,' ) = ',Poly RETURN END PROGRAM CHEBYS 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 4.6 (Chebyshev Approximation). C Section 4.5, Chebyshev Polynomials, Page 246 C PARAMETER(MaxN=15) INTEGER J,K,N REAL C,D,Sum,X,Y,Z CHARACTER ANS*1 DIMENSION T(0:MaxN,0:MaxN) DIMENSION C(0:MaxN),P(0:MaxN),X(0:MaxN),Y(0:MaxN) EXTERNAL CH,F 10 CONTINUE CALL MESSAGE(N) CALL CHEBY(C,X,Y,N) CALL CREATEP(T,N) CALL MAKECHE(T,C,P,N) CALL PRTFUN CALL PRINTCHE(C,N) CALL PRINTPOL(P,N) WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,'(A)') ANS 20 CONTINUE CALL RESULT(X,C,N) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE P(T) AGAIN ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE ANOTHER POLYNOMIAL ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END REAL FUNCTION F(X) F=EXP(X) RETURN END SUBROUTINE PRINTFUN WRITE(9,*) 'EXP(X)' RETURN END SUBROUTINE MAKECHE(T,C,P,N) PARAMETER(MaxN=15) INTEGER J,K,N REAL C,D,P,Sum,T DIMENSION T(0:MaxN,0:MaxN) DIMENSION C(0:MaxN),P(0:MaxN) DO J=0,N Sum=0.0 DO K=0,N Sum=Sum+C(K)*T(K,J) ENDDO P(J)=Sum ENDDO RETURN END SUBROUTINE CHEBY(C,X,Y,N) PARAMETER(MaxN=15) INTEGER J,K,N REAL C,D,Sum,X,Y,Z DIMENSION C(0:MaxN),X(0:MaxN),Y(0:MaxN) EXTERNAL F D=3.1415926535/(2.0*N+2.0) DO K=0,N X(K)=COS((2.0*K+1.0)*D) Y(K)=F(X(K)) C(K)=0.0 ENDDO DO K=0,N Z=(2.0*K+1.0)*D DO J=0,N C(J)=C(J)+Y(K)*COS(J*Z) ENDDO ENDDO C(0)=C(0)/(N+1.0) DO J=1,N C(J)=2.0*C(J)/(N+1.0) ENDDO RETURN END FUNCTION CH(X,C,N) PARAMETER(MaxN=15) INTEGER J,N REAL C,Sum,T,X DIMENSION C(0:MaxN),T(0:MaxN) T(0)=1.0 T(1)=X IF (N.GT.1) THEN DO J=1,N-1 T(J+1)=2.0*X*T(J)-T(J-1) ENDDO ENDIF Sum=0.0 DO J=0,N Sum=Sum+C(J)*T(J) ENDDO CH=Sum RETURN END SUBROUTINE CREATEP(T,N) PARAMETER(MaxN=15) INTEGER J,K,N REAL T DIMENSION T(0:MaxN,0:MaxN) DO K=0,N DO J=0,N T(K,J)=0.0 ENDDO ENDDO T(0,0)=1.0 T(1,0)=0.0 T(1,1)=1.0 DO K=2,N DO J=1,K T(K,J)=2.0*T(K-1,J-1)-T(K-2,J) ENDDO T(K,0)=-T(K-2,0) ENDDO RETURN END SUBROUTINE XMAKECHE(T,C,P,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL C,D,P,Sum,T DIMENSION T(0:MaxN,0:MaxN) DIMENSION C(0:MaxN),P(0:MaxN) DO 20 J=0,N Sum=0.0 DO 10 K=0,N Sum=Sum+C(K)*T(K,J) 10 CONTINUE P(J)=Sum 20 CONTINUE RETURN END SUBROUTINE XCHEBY(C,X,Y,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL C,D,Sum,X,Y,Z DIMENSION C(0:MaxN),X(0:MaxN),Y(0:MaxN) EXTERNAL F D=3.1415926535/(2.0*N+2.0) DO 10 K=0,N X(K)=COS((2.0*K+1.0)*D) Y(K)=F(X(K)) C(K)=0.0 10 CONTINUE DO 30 K=0,N Z=(2.0*K+1.0)*D DO 20 J=0,N C(J)=C(J)+Y(K)*COS(J*Z) 20 CONTINUE 30 CONTINUE C(0)=C(0)/(N+1.0) DO 40 J=1,N C(J)=2.0*C(J)/(N+1.0) 40 CONTINUE RETURN END FUNCTION XCH(X,C,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,N REAL C,Sum,T,X DIMENSION C(0:MaxN),T(0:MaxN) T(0)=1.0 T(1)=X IF (N.GT.1) THEN DO 10 J=1,N-1 T(J+1)=2.0*X*T(J)-T(J-1) 10 CONTINUE ENDIF Sum=0.0 DO 20 J=0,N Sum=Sum+C(J)*T(J) 20 CONTINUE CH=Sum RETURN END SUBROUTINE XCREATEP(T,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=15) INTEGER J,K,N REAL T DIMENSION T(0:MaxN,0:MaxN) DO 20 K=0,N DO 10 J=0,N T(K,J)=0.0 10 CONTINUE 20 CONTINUE T(0,0)=1.0 T(1,0)=0.0 T(1,1)=1.0 DO 40 K=2,N DO 30 J=1,K T(K,J)=2.0*T(K-1,J-1)-T(K-2,J) 30 CONTINUE T(K,0)=-T(K-2,0) 40 CONTINUE RETURN END SUBROUTINE MESSAGE(N) PARAMETER(MaxN=15) INTEGER N WRITE(9,*)' ' WRITE(9,*)' Construction of the Chebyshev polynomial for +the function F(X).' WRITE(9,*)' ' WRITE(9,*)' P(X) = C T (x) + C T (x) + C T (x) +...+ C T (x)' WRITE(9,*)' 0 0 1 1 2 2 N N ' WRITE(9,*)' ' WRITE(9,*)' over the interval [-1,1]. The coefficients {C } ' WRITE(9,*)' j ' WRITE(9,*)' are computed using the formulas:' WRITE(9,*)' ' WRITE(9,*)' N N' WRITE(9,*)' 1 2 2k+1 +' WRITE(9,*)' C = ---Sum f(x ) , C = ---Sum f(x )COS(J*Pi----) + for j=1,2,...,N' WRITE(9,*)' 0 N+1 k j N+1 k 2N+1' WRITE(9,*)' k=0 k=0 ' WRITE(9,*)' ' WRITE(9,*)' where x = COS(Pi*(2k+1)/(2N+2))' WRITE(9,*)' k' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ENTER the maximum number of terms N = ' N=MaxN READ(9,*) N IF (N.LT.1) THEN N=1 ENDIF IF (N.GT.MaxN) THEN N=MaxN ENDIF WRITE(9,*)' ' RETURN END SUBROUTINE PRINTPOL(A,N) PARAMETER(MaxN=15) INTEGER I,N REAL A DIMENSION A(0:MaxN) WRITE(9,*)' ' WRITE(9,*)'The Chebyshev polynomial of degree N = ',N,' is:' IF (N.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'P(X) = A X + A' WRITE(9,*)' 1 0' WRITE(9,*)' ' ELSEIF (N.EQ.2) THEN WRITE(9,*)' 2' WRITE(9,*)'P(X) = A X + A X + A' WRITE(9,*)' 2 1 0' WRITE(9,*)' ' ELSEIF (N.EQ.3) THEN WRITE(9,*)' 3 2' WRITE(9,*)'P(X) = A X + A X + A X + A ' WRITE(9,*)' 3 2 1 0' WRITE(9,*)' ' ELSE WRITE(9,*)' ',N,' 2' WRITE(9,*)'P(X) = A X +...+ A X + A X + A ' WRITE(9,*) N,' 2 1 0' WRITE(9,*)' ' ENDIF DO 10 K=N,0,-1 WRITE(9,*)'A(',K,') =',A(K) 10 CONTINUE RETURN END SUBROUTINE PRINTCHE(C,N) PARAMETER(MaxN=15) INTEGER I,K,N REAL C DIMENSION C(0:MaxN) WRITE(9,*)' ' WRITE(9,*)'The Chebyshev series with N = ',N,' terms is' IF (N.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'T(X) = C T (X) + C T (X)' WRITE(9,*)' 0 0 1 1' WRITE(9,*)' ' ELSEIF (N.EQ.2) THEN WRITE(9,*)' ' WRITE(9,*)'T(X) = C T (X) + C T (X) + C T (X)' WRITE(9,*)' 0 0 1 1 2 2' WRITE(9,*)' ' ELSEIF (N.EQ.3) THEN WRITE(9,*)' ' WRITE(9,*)'T(X) = C T (X) + C T (X) + C T (X) + C T (X)' WRITE(9,*)' 0 0 1 1 2 2 3 3' WRITE(9,*)' ' ELSE WRITE(9,*)' ' WRITE(9,*)'T(X) = C T (X) + C T (X) + ... + C T (X)' WRITE(9,*)' 0 0 1 1 N N' WRITE(9,*)' ' ENDIF DO 10 K=0,N WRITE(9,*)'C(',K,') =',C(K) 10 CONTINUE RETURN END SUBROUTINE PRTFUN WRITE(9,*) WRITE(9,*) WRITE(9,*)'You chose to construct the Chebyshev polynomial for ' WRITE(9,*) WRITE(9,*)' F(X) = EXP(X) over the interval [-1,1].' WRITE(9,*) RETURN END SUBROUTINE RESULT(X,C,N) INTEGER N REAL C,X,Z DIMENSION C(0:MaxN),X(0:MaxN) WRITE(9,*) WRITE(9,*) WRITE(9,*)'The results for the approximation of F(x) + over [-1,1] are:' WRITE(9,*) WRITE(9,*)' x F(x) P(x) + F(x) - P(x).' WRITE(9,*) DO 10 K=0,20 Z=-1 + K*0.1 WRITE(9,*) Z,' ',F(Z),' ',CH(Z,C,N), +' ',F(Z)-CH(Z,C,N) 10 CONTINUE RETURN END