skip to main |
skip to sidebar
PROGRAM PCMET
OPEN(1,'INPUT.DAT')
OPEN(2,'OUTPUT.DAT')
READ(1,*)A,B,N,Y0
CALL PCM(A,B,N,Y0)
END PROGRAM
SUBROUTINE PCM(A,B,N,Y0)
F(T,Y)=Y-T*T+1
G(T)=(T+1)**2-0.5*EXP(T)
REAL::T(0:N),W(0:N),K1,K2,K3,K4
WRITE(2,*)'ADAMS 4TH ORDER PC METHOD:'
WRITE(2,10)
10 FORMAT(3X,'T',9X,'W',8X,'EXACT',5X,'ERROR')
H=(B-A)/N
T(0)=A
W(0)=Y0
WRITE(2,11)T(0),W(0),G(T(0)),ABS(W(0)-G(T(0)))
11 FORMAT(2X,F4.2,2X,F9.6,2X,F9.6,2X,F9.6)
DO I=0,2
K1=H*F(T(I),W(I))
K2=H*F(T(I)+0.5*H,W(I)+0.5*K1)
K3=H*F(T(I)+0.5*H,W(I)+0.5*K2)
K4=H*F(T(I)+H,W(I)+K3)
W(I+1)=W(I)+(K1+2*K2+2*K3+K4)/6
T(I+1)=A+(I+1)*H
WRITE(2,11)T(I+1),W(I+1),G(T(I+1)),ABS(W(I+1)-G(T(I+1)))
END DO
DO I=3,N-1
T(I+1)=A+(I+1)*H
W(I+1)=W(I)+H*(55*F(T(I),W(I))-59*F(T(I-1),W(I-1))+37*F(T(I-2),W(
1I-2))-9*F(T(I-3),W(I-3)))/24
W(I+1)=W(I)+H*(9*F(T(I+1),W(I+1))+19*F(T(I),W(I))-5*F(T(I-1),W(I-1
1))+F(T(I-2),W(I-2)))/24
WRITE(2,11)T(I+1),W(I+1),G(T(I+1)),ABS(W(I+1)-G(T(I+1)))
END DO
END SUBROUTINE
»» read more
PROGRAM AMMET
OPEN(1,'INPUT.DAT')
OPEN(2,'OUTPUT.DAT')
READ(1,*)A,B,N,Y0
CALL AMM(A,B,N,Y0)
END PROGRAM
SUBROUTINE AMM(A,B,N,Y0)
F(T,Y)=Y-T*T+1
G(T)=(T+1)**2-0.5*EXP(T)
REAL::T(0:N),Y(0:N),W(0:N),K1,K2,K3,K4
WRITE(2,*)'AM IMPLICIT METHOD:'
WRITE(2,10)
10 FORMAT(3X,'T',9X,'W',8X,'EXACT',5X,'ERROR')
H=(B-A)/N
T(0)=A
W(0)=Y0
Y(0)=Y0
WRITE(2,11)T(0),Y(0),G(T(0)),ABS(Y(0)-G(T(0)))
11 FORMAT(2X,F4.2,2X,F9.6,2X,F9.6,2X,F9.6)
DO I=0,N-1
K1=H*F(T(I),W(I))
K2=H*F(T(I)+0.5*H,W(I)+0.5*K1)
K3=H*F(T(I)+0.5*H,W(I)+0.5*K2)
K4=H*F(T(I)+H,W(I)+K3)
W(I+1)=W(I)+(K1+2*K2+2*K3+K4)/6
T(I+1)=A+(I+1)*H
IF(I<3) THEN
Y(I+1)=W(I+1)
ELSE
Y(I+1)=Y(I)+H*(251*F(T(I+1),W(I+1))+646*F(T(I),W(I))-264*F(T(I-1),W(I-1))+106*F(T(I-2),W(I-2))-19*F(T(I-3),W(I-3)))/720
END IF
WRITE(2,11)T(I+1),Y(I+1),G(T(I+1)),ABS(Y(I+1)-G(T(I+1)))
END DO
END SUBROUTINE
»» read more
PROGRAM ABMET
OPEN(1,'INPUT.DAT')
OPEN(2,'OUTPUT.DAT')
READ(1,*)A,B,N,Y0
CALL ABM(A,B,N,Y0)
END PROGRAM
SUBROUTINE ABM(A,B,N,Y0)
F(T,Y)=Y-T*T+1
G(T)=(T+1)**2-0.5*EXP(T)
REAL::T(0:N),W(0:N),K1,K2,K3,K4
WRITE(2,*)'AB EXPLICIT METHOD:'
WRITE(2,10)
10 FORMAT(3X,'T',9X,'W',8X,'EXACT',5X,'ERROR')
H=(B-A)/N
T(0)=A
W(0)=Y0
WRITE(2,11)T(0),W(0),G(T(0)),ABS(W(0)-G(T(0)))
11 FORMAT(2X,F4.2,2X,F9.6,2X,F9.6,2X,F9.6)
DO I=0,2
K1=H*F(T(I),W(I))
K2=H*F(T(I)+0.5*H,W(I)+0.5*K1)
K3=H*F(T(I)+0.5*H,W(I)+0.5*K2)
K4=H*F(T(I)+H,W(I)+K3)
W(I+1)=W(I)+(K1+2*K2+2*K3+K4)/6
T(I+1)=A+(I+1)*H
WRITE(2,11)T(I+1),W(I+1),G(T(I+1)),ABS(W(I+1)-G(T(I+1)))
END DO
DO I=3,N-1
T(I+1)=A+(I+1)*H
W(I+1)=W(I)+H*(55*F(T(I),W(I))-59*F(T(I-1),W(I-1))+37*F(T(I-2),W(I-2))-9*F(T(I-3),W(I-3)))/24
WRITE(2,11)T(I+1),W(I+1),G(T(I+1)),ABS(W(I+1)-G(T(I+1)))
END DO
END SUBROUTINE
Input is same as RK and Euler Methods as the function, that is, initial value problem is the same. But A-B method produces more accurate approximations.
»» read more
PROGRAM RKMETS
OPEN(1,'INPUT.DAT')
OPEN(2,'OUTPUT.DAT')
READ(1,*)A,B,N,Y0
CALL RK2(A,B,N,Y0)
CALL RK4(A,B,N,Y0)
END PROGRAM
SUBROUTINE RK2(A,B,N,Y0)
F(T,Y)=Y-T*T+1
G(T)=(T+1)**2-0.5*EXP(T)
REAL::K1,K2
WRITE(2,*)'R-K METHOD OF ORDER 2:'
WRITE(2,10)
10 FORMAT(4X,"T"98X,"W",8X,'EXACT',5X,'ERROR')
H=(B-A)/N
T=A
Y=Y0
WRITE(2,11)T,Y,G(T),ABS(Y-G(T))
11 FORMAT(2X,F4.2,2X,F9.6,2X,F9.6,2X,F9.6)
DO I=1,N
K1=H*F(T,Y)
K2=H*F(T+H,Y+K1)
Y=Y+(K1+K2)/2
T=A+I*H
WRITE(2,11)T,Y,G(T),ABS(Y-G(T))
END DO
END SUBROUTINE
SUBROUTINE RK4(A,B,N,Y0)
F(T,Y)=Y-T*T+1
G(T)=(T+1)**2-0.5*EXP(T)
REAL::K1,K2,K3,K4
WRITE(2,*)'R-K METHOD OF ORDER 4:'
WRITE(2,10)
10 FORMAT(4X,"T"98X,"W",8X,'EXACT',5X,'ERROR')
H=(B-A)/N
T=A
Y=Y0
WRITE(2,11)T,Y,G(T),ABS(Y-G(T))
11 FORMAT(2X,F4.2,2X,F9.6,2X,F9.6,2X,F9.6)
DO I=1,N
K1=H*F(T,Y)
K2=H*F(T+0.5*H,Y+0.5*K1)
K3=H*F(T+0.5*H,Y+0.5*K2)
K4=H*F(T+H,Y+K3)
Y=Y+(K1+2*K2+2*K3+K4)/6
T=A+I*H
WRITE(2,11)T,Y,G(T),ABS(Y-G(T))
END DO
END SUBROUTINE
»» read more
PROGRAM ROMINT
DIMENSION R(20,20)
OPEN(1,'IA10.DAT')
OPEN(2,'OA10.DAT')
READ(1,*)A,B,N
CALL ROMBERG(A,B,N,R)
END PROGRAM
SUBROUTINE ROMBERG(A,B,N,R)
F(X)=SIN(X)
REAL::R(20,20)
WRITE(2,*)'ROMBERG INTEGRATION'
H=B-A
R(1,1)=(H/2.0)*(F(A)+F(B))
DO I=2,N
S=0
DO K=1,2**(I-2)
S=S+F(A+((K-0.5)*H))
END DO
R(I,1)=0.5*(R(I-1,1)+(H*S))
H=H/2.0
END DO
DO J=2,N
DO I=J,N
R(I,J)=R(I,J-1)+((R(I,J-1)-R(I-1,J-1))/(4**(J-1))-1)
END DO
END DO
DO I=1,N
WRITE(2,20)(R(I,J),J=1,I)
20 FORMAT(6(2X,F12.6))
END DO
END SUBROUTINE
»» read more
|