PROGRAM FINITED 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 10.1 (Finite-Difference Solution for the Wave Equation). C Section 10.1, Hyperbolic Equations, Page 507 C PARAMETER(Pi = 3.14159265) PARAMETER(MaxN = 11,MaxM = 11) CHARACTER Ans*1,Ach*1,Mess*80 REAL V DIMENSION V(1:MaxN,1:MaxM) INTEGER FunType, Inum, M, Meth, N REAL A, B, C, Rnum, Y0 EXTERNAL F,G,Fi,Gi Meth = 1 FunType = 1 A = 0 B = 1 Y0 = 0 M = 1 CALL MESSAGE(Meth) CALL INPUT(FunType) CALL EPOINTS(A,B,C,N,M,FunType) CALL FDIFF(F,G,Fi,Gi,A,B,C,N,M,U) CALL RESULT(FunType,U,N,M) END REAL FUNCTION F(X) REAL X F = SIN(3.14159265*X) + SIN(2*3.14159265*X) RETURN END REAL FUNCTION G(X) REAL X G = 0 RETURN END SUBROUTINE PRTFUN(FunType) INTEGER FunType IF (FunType.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)' The boundary functions are:' WRITE(9,*)' ' WRITE(9,*)' u(x,0) = f(x) = SIN(Pi*X) + SIN(2*Pi*X)' WRITE(9,*)' ' WRITE(9,*)' u (x,0) = g(x) = 0' WRITE(9,*)' t' ENDIF RETURN END REAL FUNCTION Fi(I,H) REAL H INTEGER I EXTERNAL F Fi = F(H*(I-1)) RETURN END REAL FUNCTION Gi(I,H) REAL H INTEGER I EXTERNAL G Gi = G(H*(I-1)) RETURN END SUBROUTINE FDIFF(F,G,Fi,Gi,A,B,C,N,M,U) PARAMETER(MaxN = 11,MaxM = 11) INTEGER N,M REAL A,B,C REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER I,J REAL H,K,R,R2,R22,S1,S2 EXTERNAL F,G,Fi,Gi H = A/(N-1) K = B/(M-1) R = C*K/H R2 = R*R R22 = R*R/2 S1 = 1 - R*R S2 = 2 - 2*R*R DO J=1,M U(1,J) = 0 U(N,J) = 0 ENDDO DO I=2,(N-1) U(I,1) = Fi(I,H) U(I,2) = S1*Fi(I,H) + K*Gi(I,H) + R22*(Fi(I+1,H) + Fi(I-1,H)) ENDDO DO J=3,M DO I=2,(N-1) U(I,J) = S2*U(I,J-1) + R2*(U(I-1,J-1) + U(I+1,J-1)) - U(I,J-2) ENDDO ENDDO RETURN END SUBROUTINE MESSAGE(Meth) INTEGER K,Meth WRITE(9,*)' SOLUTION OF HYPERBOLIC EQUATIONS' WRITE(9,*)' ' Meth = 1 RETURN END SUBROUTINE INPUT(FunType) INTEGER K,FunType CHARACTER Ans*1 WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'The finite difference method is used' WRITE(9,*)'to solve the wave equation' WRITE(9,*)' ' WRITE(9,*)' 2 ' WRITE(9,*)' u (x,t) = c u (x,t)' WRITE(9,*)' tt xx' WRITE(9,*)' ' WRITE(9,*)'with u(0,t) = 0 and u(a,t) = 0 for 0 ² t ² B.' WRITE(9,*)' ' WRITE(9,*)'and' WRITE(9,*)' ' WRITE(9,*)'u(x,0) = f(x) and u (x,0) = g(x) for 0 < x < A.' WRITE(9,*)' t' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'A numerical approximation is' WRITE(9,*)'computed over the rectangle' WRITE(9,*)' ' WRITE(9,*)' 0 ² x ² A.' WRITE(9,*)' 0 ² t ² B.' WRITE(9,*)' ' WRITE(9,*)'You must supply the endpoints for the intervals.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Press the key. ' READ(9,*) Ans WRITE(9,*)' ' FunType = 1 RETURN END SUBROUTINE EPOINTS(A, B, C, N, M,FunType) INTEGER I,N,M,FunType REAL A, B, C, Valu CHARACTER Resp*1 CHARACTER Mess*80 WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunType) WRITE(9,*)' ' WRITE(9,*)' ' Mess = 'For the interval [0,A], ENTER the endpoint A = ' WRITE(9,*) Mess READ(9,*) A WRITE(9,*)' ' Mess = 'For the interval [0,B], ENTER the endpoint B = ' WRITE(9,*) Mess READ(9,*) B WRITE(9,*)' ' Mess = ' ENTER the constant C = ' WRITE(9,*) Mess READ(9,*) C WRITE(9,*)' ' Mess = ' ENTER the number of steps N = ' WRITE(9,*) Mess READ(9,*) N IF (N.LT.2) THEN N = 2 ENDIF IF (N.GT.25) THEN N = 25 ENDIF WRITE(9,*)' ' Mess = ' ENTER the number of steps M = ' WRITE(9,*) Mess READ(9,*) M IF (M.LT.2) THEN M = 2 ENDIF IF (M.GT.100) THEN M = 25 ENDIF RETURN END SUBROUTINE RESULT(FunType,U,N,M) PARAMETER(MaxN = 11,MaxM = 11) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER FunType,I,J,N,M WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunType) WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' u(x ,t ) ..... u(x ,t )' WRITE(9,*)' 2 j N-1 j' WRITE(9,*)' -------------------------------------------------' WRITE(9,*)' ' DO J=1,M WRITE(9,999) (U(I,J), I=2,(N-1)) ENDDO WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,*) Ans WRITE(9,*)' ' 999 FORMAT(1X,9F8.4) RETURN END PROGRAM HEATEQN 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 10.2 (Forward-Difference Method for the Heat Equation). C Section 10.2, Parabolic Equations, Page 516 C C Algorithm 10.3 (Crank-Nicholson Method for the Heat Equation). C Section 10.2, Parabolic Equations, Page 517 C C PRENTICE HALL, INC. © 1992 PARAMETER(Pi = 3.1415926535) PARAMETER(MaxN = 11,MaxM = 11) CHARACTER Ans*1,Ach*1,Mess*80 REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER FunType, Inum, M, Meth, N REAL A, B, C, C1, C2, Rnum, Y0 EXTERNAL F,G1,G2,Fi,G1i,G2i Meth = 1 FunType = 1 A = 0 B = 1 Y0 = 0 M = 1 CALL MESSAGE(Meth) CALL INPUT(FunType,Meth) CALL EPOINTS(A,B,C,N,M,FunType) IF (Meth.EQ.1) THEN CALL ForDiff(F,G1,G2,A,B,C,N,M,U) ENDIF IF (Meth.EQ.2) THEN CALL CraNich(F,G1,G2,A,B,C,N,M,U) ENDIF CALL RESULT(FunType,U,N,M) END REAL FUNCTION F(X) REAL X F = SIN(3.1415926535 * X) + SIN(3 * 3.1415926535 * X) RETURN END REAL FUNCTION G1(T) REAL T G1 = 0 RETURN END REAL FUNCTION G2(T) REAL T G2 = 0 RETURN END SUBROUTINE PRTFUN(FunType) INTEGER FunType WRITE(9,*)' ' WRITE(9,*)' The boundary functions are:' WRITE(9,*)' ' WRITE(9,*)' u(x,0) = f(x) = SIN(Pi*X) + SIN(3*Pi*X)' WRITE(9,*)' ' WRITE(9,*)' u(0,t) = g1(t) = 0' WRITE(9,*)' ' WRITE(9,*)' u(a,t) = g2(t) = 0' RETURN END REAL FUNCTION Fi(I,H) REAL H INTEGER I Fi = F(H * (I - 1)) RETURN END REAL FUNCTION G1i(I,K) REAL K INTEGER I G1i = G1(K * (I - 1)) RETURN END REAL FUNCTION G2i(I,K) REAL K INTEGER I G2i = G2(K * (I - 1)) RETURN END SUBROUTINE ForDiff(F,G1,G2,A,B,C,N,M,U) PARAMETER(MaxN = 11,MaxM = 11) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER I,J,N,M REAL A,B,C,H,K,R,S EXTERNAL F, G1, G2, Fi, G1i, G2i H = A/(N-1) K = B/(M-1) R = C*C*K/H/H S = 1 - 2 * R DO J=1,M U(1,J) = G1i(J,K) U(N,J) = G2i(J,K) ENDDO DO I=2,(N-1) U(I,1) = Fi(I,H) ENDDO DO J=2,M DO I=2,(N-1) U(I,J) = S*U(I,J-1) + R*(U(I-1,J-1) + U(I+1,J-1)) ENDDO ENDDO RETURN END SUBROUTINE TriSys(Va,Vb,Vc,Vd,X0,N) PARAMETER(MaxN = 11) REAL X0,Va,Vb,Vc,Vd DIMENSION X0(1:MaxN),Va(1:MaxN),Vb(1:MaxN),Vc(1:MaxN),Vd(1:MaxN) REAL A0,B0,C0,D0 DIMENSION A0(1:MaxN),B0(1:MaxN),C0(1:MaxN),D0(1:MaxN) INTEGER K,N REAL T DO K=1,N A0(K) = Va(K) B0(K) = Vb(K) C0(K) = Vc(K) D0(K) = Vd(K) ENDDO DO K=2,N T = A0(K-1)/D0(K-1) D0(K) = D0(K) - T*C0(K-1) B0(K) = B0(K) - T*B0(K-1) ENDDO X0(N) = B0(N) / D0(N) DO K=(N-1),1,-1 X0(K) = (B0(K) - C0(K) * X0(K+1))/D0(K) ENDDO RETURN END SUBROUTINE CraNich(F,G1,G2,A,B,C,N,M,U) PARAMETER(MaxN = 11,MaxM = 11) REAL X,X0,Va,Vb,Vc,Vd DIMENSION X(1:MaxN),X0(1:MaxN),Va(1:MaxN) DIMENSION Vb(1:MaxN),Vc(1:MaxN),Vd(1:MaxN) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER I,J,N,M REAL A,B,C,H,K,R,S1,S2 EXTERNAL F,G1,G2 H = A/(N-1) K = B/(M-1) R = C*C*K/H/H S1 = 2 + 2/R S2 = 2/R - 2 DO J=1,M U(1,J) = G1i(J,K) U(N,J) = G2i(J,K) ENDDO DO I=2,(N-1) U(I,1) = Fi(I,H) ENDDO DO I=1,N Vd(I) = S1 ENDDO Vd(1) = 1 Vd(N) = 1 DO I=1,(N-1) Va(I) = -1 Vc(I) = -1 ENDDO Va(N - 1) = 0 Vc(1) = 0 Vb(1) = G1i(1,K) Vb(N) = G2i(1,K) DO J=2,M DO I=2,(N-1) Vb(I) = U(I-1,J-1) + U(I+1,J-1) + S2*U(I,J-1) ENDDO CALL TriSys(Va,Vb,Vc,Vd,X,N) DO I=1,N U(I,J) = X(I) ENDDO ENDDO RETURN END SUBROUTINE MESSAGE(Meth) INTEGER K,Meth WRITE(9,*)' SOLUTION OF PARABOLIC EQUATIONS' WRITE(9,*)' ' Meth = 1 RETURN END SUBROUTINE INPUT(FunType,Meth) INTEGER K,FunType,Meth CHARACTER Ans*1,Mess*80 WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Solution of the heat equation' WRITE(9,*)' ' WRITE(9,*)' 2 ' WRITE(9,*)' u (x,t) = c u (x,t)' WRITE(9,*)' t xx' WRITE(9,*)' ' WRITE(9,*)'and u(x,0) = f(x) for 0 < x < A.' WRITE(9,*)' ' WRITE(9,*)'and' WRITE(9,*)' ' WRITE(9,*)'u(0,t) = g1(t) and u(a,t) = g2(t) for 0 ² t ² B.' WRITE(9,*)' ' FunType = 1 WRITE(9,*)'A numerical approximation is' WRITE(9,*)'computed over the rectangle' WRITE(9,*)' ' WRITE(9,*)' 0 ² x ² A.' WRITE(9,*)' 0 ² t ² B.' WRITE(9,*)' ' WRITE(9,*)'You must supply the endpoints for the intervals.' WRITE(9,*)' ' WRITE(9,*)' Choose the method of computation:' WRITE(9,*)' ' WRITE(9,*)' < 1 > The forward difference method.' WRITE(9,*)' ' WRITE(9,*)' < 2 > The Crank-Nicholson method.' WRITE(9,*)' ' Mess = ' SELECT your method < 1 or 2 > ? ' WRITE(9,*) Mess READ(9,*) Meth IF (Meth.LT.1) THEN Meth = 1 ENDIF IF (Meth.GT.2) THEN Meth = 2 ENDIF WRITE(9,*)' ' RETURN END SUBROUTINE EPOINTS(A,B,C,N,M,FunType) INTEGER I,N,M REAL A,B,C,Valu CHARACTER Resp*1 CHARACTER Ans*1,Ach*1,Mess*80 CALL PRTFUN(FunType) WRITE(9,*)' ' Mess = 'For the interval [0,A], ENTER the endpoint A = ' WRITE(9,*) Mess READ(9,*) A WRITE(9,*)' ' Mess = 'For the interval [0,B], ENTER the endpoint B = ' WRITE(9,*) Mess READ(9,*) B WRITE(9,*)' ' Mess = ' ENTER the constant C = ' WRITE(9,*) Mess READ(9,*) C WRITE(9,*)' ' Mess = ' ENTER the number of steps N = ' WRITE(9,*) Mess READ(9,*) N IF (N.LT.2) THEN N = 2 ENDIF IF (N.GT.25) THEN N = 25 ENDIF WRITE(9,*)' ' Mess = ' ENTER the number of steps M = ' WRITE(9,*) Mess READ(9,*) M IF (M.LT.2) THEN M = 2 ENDIF IF (M.GT.100) THEN M = 25 ENDIF RETURN END SUBROUTINE RESULT(FunType,U,N,M) PARAMETER(MaxN = 11,MaxM = 11) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER FunType,I,J,N,M WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunType) WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' u(x ,t ) ..... u(x ,t )' WRITE(9,*)' 2 j N-1 j' WRITE(9,*)' ---------------------------------------------------' WRITE(9,*)' ' DO J=1,M WRITE(9,999) (U(I,J), I=2,(N-1)) ENDDO WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,*) Ans WRITE(9,*)' ' 999 FORMAT(1X,9F8.4) RETURN END PROGRAM LAPLACE 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 10.4 (Dirichlet Method for Laplace's Equation). C Section, 10.3 Elliptic Equations, Page 531 C PARAMETER(Pi = 3.1415926535) PARAMETER(MaxN = 11,MaxM = 11) CHARACTER Ans*1,Ach*1,Mess*80 REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER FunType, M, Meth, N REAL A, B, C, C1, C2, H, Rnum, Y0 EXTERNAL F1,F2,F3,F4,F1i,F2i,F3i,F4i Meth = 1 FunType = 1 A = 0 B = 1 Y0 = 0 M = 1 CALL MESSAGE(Meth) CALL INPUT(FunType, Meth) CALL EPOINTS(A,B,H,N,M) CALL Dirich(F1,F2,F3,F4,A,B,H,N,M,U) CALL RESULTS(FunType,U,N,M) END REAL FUNCTION F1(X) REAL X F1 = 20 RETURN END REAL FUNCTION F2(X) REAL X F2 = 180 RETURN END REAL FUNCTION F3(X) REAL X F3 = 80 RETURN END REAL FUNCTION F4(X) REAL X F4 = 0 RETURN END SUBROUTINE PRTFUN(FunType) INTEGER FunType WRITE(9,*)' ' WRITE(9,*)' The boundary functions are:' WRITE(9,*)' ' WRITE(9,*)' u(x,0) = f1(x) = 20' WRITE(9,*)' ' WRITE(9,*)' u(x,b) = f2(x) = 180' WRITE(9,*)' ' WRITE(9,*)' u(0,y) = f3(x) = 80' WRITE(9,*)' ' WRITE(9,*)' u(a,y) = f4(x) = 0' RETURN END REAL FUNCTION F1i(I,H) INTEGER I REAL H F1i = F1(H*(I-1)) RETURN END REAL FUNCTION F2i(I,H) INTEGER I REAL H F2i = F2(H*(I-1)) RETURN END REAL FUNCTION F3i(I,H) INTEGER I REAL H F3i = F3(H*(I-1)) RETURN END REAL FUNCTION F4i(I,H) INTEGER I REAL H F4i = F4(H*(I-1)) RETURN END SUBROUTINE Dirich(F1,F2,F3,F4,A,B,H,N,M,U) PARAMETER(Pi = 3.1415926535) PARAMETER(MaxN = 11,MaxM = 11) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER Count,I,J,N,M REAL A,B,H,Ave,K,R,S,W,Relax,Tol EXTERNAL F1,F2,F3,F4,F1i,F2i,F3i,F4i Ave = (A*(f1(0)+f2(0)) + B*(f3(0)+f4(0)))/(2*A+2*B) DO I=2,(N-1) DO J=2,(M-1) U(I,J) = Ave ENDDO ENDDO DO J=1,M U(1,J) = F3i(J) U(N,J) = F4i(J) ENDDO DO I=1,N U(I,1) = F1i(I) U(I,M) = F2i(I) ENDDO U(1,1) = (U(1,2) + U(2,1))/2 U(1,M) = (U(1,M-1) + U(2,M))/2 U(N,1) = (U(N-1,1) + U(N,2))/2 U(N,M) = (U(N-1,M) + U(N,M-1))/2 W = 4/(2+SQRT(4-(COS(Pi/(N-1))+COS(Pi/(M-1)))**2)) Tol = 1 Count = 0 WHILE ((Tol.GT.0.001).AND.(Count.LE.25)) Tol = 0 DO J=2,(M-1) DO I=2,(N-1) Relax = W*(U(I,J+1) + U(I,J-1) + U(I+1,J) + & U(I-1,J) - 4.0*U(I,J))/4.0 U(I,J) = U(I,J) + Relax IF (Tol.LE.ABS(Relax)) THEN Tol = ABS(Relax) ENDIF ENDDO ENDDO Count = Count+1 REPEAT RETURN END SUBROUTINE MESSAGE(Meth) INTEGER K,Meth WRITE(9,*)' SOLUTION OF ELLIPTIC EQUATIONS' WRITE(9,*)' ' Meth = 1 RETURN END SUBROUTINE INPUT(FunType,Meth) INTEGER K,FunType,Meth CHARACTER Ans*1 WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Solution of Laplace`s equation' WRITE(9,*)' ' WRITE(9,*)' u (x,y) + u (x,y) = 0' WRITE(9,*)' xx yy' WRITE(9,*)' ' WRITE(9,*)'with the boundary values:' WRITE(9,*)' ' WRITE(9,*)'u(x,0)=f1(x), u(x,b)=f2(x) for 0 ² x ² A,' WRITE(9,*)' ' WRITE(9,*)'u(0,y)=f3(x), u(a,y)=f4(x) for 0 ² x ² B,' WRITE(9,*)' ' WRITE(9,*)'A numerical approximation is' WRITE(9,*)'computed over the rectangle' WRITE(9,*)' ' WRITE(9,*)' 0 ² x ² A.' WRITE(9,*)' 0 ² t ² B.' WRITE(9,*)' ' WRITE(9,*)'You supply the endpoints for the intervals.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'Press the key. ' READ(9,*) Ans WRITE(9,*)' ' FunType = 1 Meth = 1 RETURN END SUBROUTINE EPOINTS(A,B,H,N,M) INTEGER I,N,M REAL A, B, H,Valu CHARACTER Resp*1,Mess*80 WRITE(9,*)' ' CALL PRTFUN(FunType) WRITE(9,*)' ' Mess = 'For the interval[0,A], ENTER the endpoint A =' WRITE(9,*) Mess READ(9,*) A WRITE(9,*)' ' Mess = 'For the interval[0,B], ENTER the endpoint B =' WRITE(9,*) Mess READ(9,*) B WRITE(9,*)' ' Mess = ' ENTER the step size H =' WRITE(9,*) Mess READ(9,*) H H = ABS(H) N = ANINT(A/H)+1 M = ANINT(B/H)+1 RETURN END SUBROUTINE RESULTS(FunType,U,N,M) PARAMETER(MaxN = 11,MaxM = 11) REAL U DIMENSION U(1:MaxN,1:MaxM) INTEGER FunType,I,J,N,M WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRTFUN(FunType) WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' u(x ,y ) ..... u(x ,y )' WRITE(9,*)' 1 j N j' WRITE(9,*)' ---------------------------------------------' WRITE(9,*)' ' DO J=M,1,-1 WRITE(9,999) (U(I,J), I=1,N) ENDDO WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,*) Ans WRITE(9,*)' ' 999 FORMAT(1X,9F8.3) RETURN END