Wednesday, March 17, 2010

Assignment 15 : ADAMS 4TH ORDER PREDICTOR-CORRECTOR METHOD

      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

Assignment 14 : Adam - Moulton 4 Step Implicit Method

      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

Assignment 13 : Adam - Bashforth 4 Step Explicit Method

      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

Assignment 12 : RK Method of Order 2 & 4

      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

Assignment 10 : Romberg Integration

      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