MODULE const INTEGER, PARAMETER :: dp = KIND(1.0d0), sp = KIND(1.0e0) REAL(KIND=dp), PARAMETER :: dpi = 3.141592653589793d0 END MODULE const !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PROGRAM eulertest USE const IMPLICIT NONE ! ! Program to implement Adams 4th order predictor-corrector method for the problem ! ! y' = -30y + cos(t), y(0) = 3 on the interval [0 4] ! INTEGER,PARAMETER :: NMAX = 2000001 INTEGER :: n,i REAL(kind=dp) :: a,b,ic,tn REAL(kind=dp) :: t(NMAX),y(NMAX),ye(NMAX),error(NMAX) REAL(kind=dp) :: norm ! Set up inputs (ic is the initial condition) a = 0 b = 4 ic = 3 WRITE(*,*) 'Input n' READ(*,*) n ! Call the main Euler routine ! n,a,b,ic are inputs ! t and y are outputs CALL ab4(n,a,b,ic,t,y) ! Compute the exact solution and error vector DO i = 1,n+1 ye(i) = (2673*exp(-30*t(i)) + sin(t(i)) + 30*cos(t(i)))/901 error(i) = y(i)-ye(i) ENDDO ! Get 2 norm of error vector tn = norm(n+1,error) WRITE(*,*) '2 norm error = ',tn END PROGRAM eulertest !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE ab4(n,a,b,ic,tvec,yvec) USE const IMPLICIT NONE ! ! Subroutine to implement 4th order Adams predictor-corrector ! method ! ic = initial condition ! tvec = vector of time slices ! yvec = approximate solution at tvec values ! INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: a,b,ic REAL(kind=dp),INTENT(OUT) :: tvec(n+1),yvec(n+1) INTEGER :: i REAL(kind=dp) :: h,fvec(n+1),k1,k2,k3,k4,temp,temp2 REAL(kind=dp) :: rhsfun ! Get t coordinates and h CALL getxcoord(n,a,b,tvec) h = tvec(2)-tvec(1) ! Set initial condition for yvec, get first function yvec(1) = ic fvec(1) = rhsfun(tvec(1),yvec(1)) ! Need to generate next 3 points using RK4 DO i = 1,3 k1 = fvec(i) k2 = rhsfun(tvec(i)+h/2,yvec(i)+h*k1/2) k3 = rhsfun(tvec(i)+h/2,yvec(i)+h*k2/2) k4 = rhsfun(tvec(i+1),yvec(i)+h*k3) yvec(i+1) = yvec(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6 fvec(i+1) = rhsfun(tvec(i+1),yvec(i+1)) ENDDO ! Remaining steps using ab4 DO i = 4,n temp = yvec(i) + h*(55*fvec(i)-59*fvec(i-1)+37*fvec(i-2)-9*fvec(i-3))/24 temp2 = rhsfun(tvec(i+1),temp) yvec(i+1) = yvec(i) + h*(9*temp2+19*fvec(i)-5*fvec(i-1)+fvec(i-2))/24 fvec(i+1) = rhsfun(tvec(i+1),yvec(i+1)) ENDDO END SUBROUTINE ab4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION rhsfun(tval,yval) USE const IMPLICIT NONE ! ! Function to evaluate the right hand function of a ! first order initial value problem. ! REAL(kind=dp),INTENT(IN) :: tval,yval REAL(kind=dp) :: rhsfun rhsfun = -30*yval + COS(tval) END FUNCTION rhsfun !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION norm(n,x) USE const IMPLICIT NONE ! ! Function to get the 2 norm of a vector ! INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: x(n) REAL(kind=dp) :: norm INTEGER :: i norm = 0 DO i = 1,n norm = norm + x(i)**2 ENDDO norm = SQRT(norm) END FUNCTION norm !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE getxcoord(n,a,b,x) ! ! Subroutine to generate a set of equally spaced ! x coordinates given the interval limits [a,b] and ! the number of subdivisions n ! USE const IMPLICIT NONE INTEGER,INTENT(IN) :: n REAL(kind=dp),INTENT(IN) :: a,b REAL(kind=dp),INTENT(OUT) :: x(n+1) INTEGER :: i REAL(kind=dp) :: h h = (b-a)/n DO i = 1,n+1 x(i) = a +(i-1)*h ENDDO RETURN END SUBROUTINE getxcoord