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 Euler's method for the problem ! ! y' = -ty, y(0) = 2 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 = 2 WRITE(*,*) 'Input n' READ(*,*) n ! Call the main Euler routine ! n,a,b,ic are inputs ! t and y are outputs CALL euler(n,a,b,ic,t,y) ! Compute the exact solution and error vector DO i = 1,n+1 ye(i) = 2*EXP(-t(i)**2/2) error(i) = y(i)-ye(i) ENDDO ! Get 2 norm of error vector tn = norm(n+1,error) DO i = n+1-10,n+1 WRITE(*,*) t(i),y(i),ye(i) ENDDO WRITE(*,*) '2 norm error = ',tn END PROGRAM eulertest !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE euler(n,a,b,ic,tvec,yvec) USE const IMPLICIT NONE ! ! Subroutine to implement Euler's 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 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 yvec(1) = ic ! Main Euler loop DO i = 1,n yvec(i+1) = yvec(i) + h*rhsfun(tvec(i),yvec(i)) ENDDO END SUBROUTINE euler !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 = -tval*yval 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