Solution to the Test Exam
!_______________________________________________________________________
! Donev Aleksandar, PHY201, a Fortran 90 program
! that implements the trapezoid rule for numerical integration
! and integrates the function f(x)=4/(1+x^2) on the interval [0,1].
! The analytical result is known to be Pi.
!_______________________________________________________________________
!_______________________________________________________________________
! This module contains single and double precision kind constants for reals
!_______________________________________________________________________
module Precision
implicit none
public
integer, parameter :: sp=kind(0.0e0), dp=kind(0.0d0) ! Single and double
integer, parameter :: wp=sp ! Or dp
end module Precision
!_______________________________________________________________________
! This module implements the trapezoid rule for numerical integration
! of a smooth function f(x) over a finite interval [a,b] with
! a division with N intervals:
!
! Integrate[f(x),x=a..b] <> h*Sum[f(a+i*h),i=1..(N-1)]+(f(a)+f(b))/2
! where h=|b-a|/N
! Modeled after the solution of Khoa Lam for worksheet 8
!_______________________________________________________________________
MODULE Trapzd
USE Precision
PUBLIC
CONTAINS
!___________________________________
! This function integrates the
! function f(x) over [a,b]
! using N intervals
! in the trapezoid rule:
!___________________________________
FUNCTION Trapezoid(f,a,b,N) RESULT(integral)
!_______________
! Input arguments and result declarations:
!_______________
INTERFACE ! Interface for function to be integrated
FUNCTION f(x) RESULT(f_x)
USE Precision
REAL(KIND=wp), INTENT(IN) :: x
REAL(KIND=wp) :: f_x
END FUNCTION f
END INTERFACE
REAL(KIND=wp), INTENT(IN) :: a,b ! The bounds
INTEGER, INTENT(IN) :: N
REAL(KIND=wp) :: integral ! The result
!_______________
! Local variables declarations:
!_______________
REAL(KIND=wp) :: h, sum ! Step size and temporary for the sum
INTEGER :: i ! Just a counter
!_______________
! Calculation:
!_______________
h=(b-a)/N
! Calculate the sum:
sum=0.0_wp ! First initialize
DO i=1, N-1
sum=sum+f(a+i*h) ! Or use x=x+h, starting at x=a
END DO
! The result is:
integral=h*sum+(h/2)*(f(a)+f(b))
END FUNCTION Trapezoid
END MODULE Trapzd
!_______________________________________________________________________
! This module contains the function f(x)=4/(1+x^2)
!_______________________________________________________________________
module Function
use Precision
implicit none
public
contains
!___________________________________
! This is the function f(x)
!___________________________________
function f(x)
real(kind=wp), intent(in) :: x
real(kind=wp) :: f
f=4.0_wp/(1.0_wp+x**2)
end function f
end module Function
!_______________________________________________________________________
!_______________________________________________________________________
! The main program calls Trapzd passing f to it to obtain the integral:
! Integrate[1/(1+x^2),x=0..1]=Pi and writes the result to the file Pi.dat
!_______________________________________________________________________
program Pi
use Precision
use Trapzd
use Function
implicit none
integer :: N ! The number of intervals
! Open the file:
open(unit=10,file="Pi.dat",action="write",status="unknown")
! Do the calculation
write(unit=*,fmt=*) "How many intervals to use?" ! Ask the user
read(unit=*,fmt=*) N ! The user enters N
! Call Trapzd and write the result to the file
write(unit=10,fmt=*) "The integral value is: ",Trapezoid(f,0.0_wp,1.0_wp,N)
write(unit=10,fmt=*) "True value of Pi is: ",4.0_wp*atan(1.0_wp)
! Close the file:
close(unit=10)
end program Pi