! !_______________________________________________________________________ ! 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 !