Implemented Several Numerical Integration Algorithms in maths/

DIRECTORY.md

@@ -1,5 +1,11 @@
 # The Algorithms - Directory of Fortran
 ## Maths
+* numerical_integration
+  * [trapezoidal_rule](/modules/maths/numerical_integration/trapezoid.f90)
+  * [simpson_rule](/modules/maths/numerical_integration/simpson.f90)
+  * [midpoint_rule](/modules/maths/numerical_integration/midpoint.f90)
+  * [monte_carlo](/modules/maths/numerical_integration/monte_carlo.f90)
+  * [gauss_legendre](/modules/maths/numerical_integration/gaussian_legendre.f90)
 * [euclid_gcd](/modules/maths/euclid_gcd.f90)
 * [factorial](/modules/maths/factorial.f90)
 * [fibonacci](/modules/maths/fibonacci.f90)

examples/maths/numerical_integration/gaussian_legendre.f90

@@ -0,0 +1,37 @@
+!> Example Program for Gaussian-Legendre Quadrature Module
+!! This program demonstrates the use of Gaussian-Legendre Quadrature Module for numerical integration.
+!!
+!! It sets the integration limits and the number of quadrature points (n), and calls the
+!! gauss_legendre_quadrature subroutine to compute the approximate value of the definite integral
+!! of the specified function.
+!!
+!! Example function: f(x) = exp(-x^2) * cos(2.0_dp * x)
+
+program example_gaussian_quadrature
+    use gaussian_legendre_quadrature
+    implicit none
+
+    real(dp) :: a, b, integral_result
+    integer :: n
+
+    ! Set the integration limits and number of quadrature points
+    a = -1.0_dp
+    b = 1.0_dp
+    n = 5       !! Number of quadrature points
+
+    ! Call Gaussian quadrature to compute the integral
+    call gauss_legendre_quadrature(integral_result, a, b, n, func)
+
+    write(*, '(A, F12.6)') "Gaussian Quadrature result: ", integral_result          !! ≈ 0.858574
+
+contains
+
+    function func(x) result(fx)
+        implicit none
+        real(dp), intent(in) :: x
+        real(dp) :: fx
+
+        fx = exp(-x**2) * cos(2.0_dp * x)       !! Example function to integrate
+    end function func
+
+end program example_gaussian_quadrature

examples/maths/numerical_integration/midpoint.f90

@@ -0,0 +1,36 @@
+!> Example Program for Midpoint Rule
+!! This program demonstrates the use of Midpoint Rule for numerical integration.
+!!
+!! It sets the integration limits and number of subintervals (panels), and calls the
+!! midpoint subroutine to compute the approximate value of the definite integral
+!! of the specified function.
+!!
+!! Example function: f(x) = exp(-x^2) * cos(2.0_dp * x)
+
+program example_midpoint
+    use midpoint_rule
+    implicit none
+    real(dp) :: a, b, integral_result
+    integer :: n
+
+    ! Set the integration limits and number of panels
+    a = -1.0_dp
+    b = 1.0_dp
+    n = 400     !! Number of subdivisions
+
+    ! Call the midpoint rule subroutine with the function passed as an argument
+    call midpoint(integral_result, a, b, n, func)
+
+    write(*, '(A, F12.6)') "Midpoint rule yields: ", integral_result  !! ≈ 0.858196
+
+contains
+
+    function func(x) result(fx)
+        implicit none
+        real(dp), intent(in) :: x
+        real(dp) :: fx
+
+        fx = exp(-x**2) * cos(2.0_dp * x)       !! Example function to integrate
+    end function func
+
+end program example_midpoint

examples/maths/numerical_integration/monte_carlo.f90

@@ -0,0 +1,37 @@
+!> Example Program for Monte Carlo Integration
+!! This program demonstrates the use of Monte Carlo module for numerical integration.
+!!
+!! It sets the integration limits and number of random samples, and calls the
+!! monte_carlo subroutine to compute the approximate value of the definite integral
+!! of the specified function.
+!!
+!! Example function: f(x) = exp(-x^2) * cos(2.0_dp * x)
+
+program example_monte_carlo
+    use monte_carlo_integration
+    implicit none
+
+    real(dp) :: a, b, integral_result, error_estimate
+    integer :: n
+
+    ! Set the integration limits and number of random samples
+    a = -1.0_dp
+    b = 1.0_dp
+    n = 1E6     !! Number of random samples
+
+    ! Call Monte Carlo integration
+    call monte_carlo(integral_result, error_estimate, a, b, n, func)
+
+    write(*, '(A, F12.6, A, F12.6)') "Monte Carlo result: ", integral_result, " +- ", error_estimate     !! ≈ 0.858421
+
+contains
+
+    function func(x) result(fx)
+        implicit none
+        real(dp), intent(in) :: x
+        real(dp) :: fx
+
+        fx = exp(-x**2) * cos(2.0_dp * x)       !! Example function to integrate
+    end function func
+
+end program example_monte_carlo

examples/maths/numerical_integration/simpson.f90

@@ -0,0 +1,37 @@
+!> Example Program for Simpson's Rule
+!! This program demonstrates the use of Simpson's rule for numerical integration.
+!!
+!! It sets the integration limits and number of panels, and calls the
+!! simpson subroutine to compute the approximate value of the definite integral
+!! of the specified function.
+!!
+!! Example function: f(x) = exp(-x^2) * cos(2.0_dp * x)
+
+program example_simpson
+    use simpson_rule
+    implicit none
+
+    real(dp) :: a, b, integral_result
+    integer :: n
+
+    ! Set the integration limits and number of panels
+    a = -1.0_dp
+    b = 1.0_dp
+    n = 100     !! Number of subdivisions (must be even)
+
+    ! Call Simpson's rule with the function passed as an argument
+    call simpson(integral_result, a, b, n, func)
+
+    write(*, '(A, F12.8)') "Simpson's rule yields: ", integral_result       !! ≈ 0.85819555
+
+contains
+
+    function func(x) result(fx)
+        implicit none
+        real(dp), intent(in) :: x
+        real(dp) :: fx
+
+        fx = exp(-x**2) * cos(2.0_dp * x)       !! Example function to integrate
+    end function func
+
+end program example_simpson

examples/maths/numerical_integration/trapezoid.f90

@@ -0,0 +1,37 @@
+!> Example Program for Trapezoidal Rule
+!! This program demonstrates the use of the Trapezoidal rule for numerical integration.
+!!
+!! It sets the integration limits and number of panels, and calls the
+!! trapezoid subroutine to compute the approximate value of the definite integral
+!! of the specified function.
+!!
+!! Example function: f(x) = exp(-x^2) * cos(2.0_dp * x)
+
+program example_tapezoid
+    use trapezoidal_rule
+    implicit none
+
+    real(dp) :: a, b, integral_result
+    integer :: n
+
+    ! Set the integration limits and number of panels
+    a = -1.0_dp
+    b = 1.0_dp
+    n = 1E6     !! Number of subdivisions
+
+    ! Call the trapezoidal rule with the function passed as an argument
+    call trapezoid(integral_result, a, b, n, func)
+
+    write(*, '(A, F12.6)') 'Trapezoidal rule yields: ', integral_result     !! ≈ 0.858195
+
+contains
+
+    function func(x) result(fx)
+        implicit none
+        real(dp), intent(in) :: x
+        real(dp) :: fx
+
+        fx = exp(-x**2) * cos(2.0_dp * x)       !! Example function to integrate
+    end function func
+
+end program example_tapezoid

modules/maths/numerical_integration/gaussian_legendre.f90

@@ -0,0 +1,100 @@
+!> Gaussian-Legendre Quadrature Module
+!!
+!! This module provides the implementation of Gaussian-Legendre Quadrature.
+!!
+!! The method approximates the definite integral of a function over a specified interval [a, b].
+!!
+!! The quadrature method works by transforming nodes and weights from the reference interval [-1, 1] to the
+!! interval [a, b] and then evaluating the function at these nodes. The integral is then approximated by summing
+!! the product of function values and corresponding weights.
+!!
+!! Contents:
+!! - `gauss_legendre_quadrature`: A subroutine to perform Gaussian-Legendre quadrature using provided nodes and weights.
+!! - `gauss_legendre_weights`: A helper subroutine to initialize the quadrature nodes and weights for different orders (n).
+!!
+!! Input:
+!! - `a`: Lower bound of integration (real(dp))
+!! - `b`: Upper bound of integration (real(dp))
+!! - `n`: Number of quadrature points (integer)
+!! - `func`: The function to integrate (interface)
+!!
+!! Output:
+!! - `integral_result`: Approximate value of the integral (real(dp))
+
+module gaussian_legendre_quadrature
+    implicit none
+    integer, parameter :: dp = kind(1.0d0)      !! Double precision parameter
+
+contains
+
+    ! General Gaussian Quadrature for definite integral
+    subroutine gauss_legendre_quadrature(integral_result, a, b, n, func)
+        implicit none
+        real(dp), intent(out) :: integral_result
+        real(dp), intent(in) :: a, b
+        integer, intent(in) :: n         !! Number of quadrature points (order of accuracy)
+
+        real(dp), dimension(n) :: t, w, x
+        real(dp), dimension(:), allocatable :: fx
+        integer :: i
+
+        ! Interface for the function
+        interface
+            real(kind(0.d0)) function func(x) result(fx)
+                real(kind(0.d0)), intent(in) :: x
+            end function func
+        end interface
+
+        ! Initialize nodes and weights for Gauss-Legendre quadrature based on n
+        call gauss_legendre_weights(t, w, n)
+
+        ! Allocate the function value array
+        allocate(fx(n))
+
+        ! Transform the nodes from the reference interval [-1, 1] to [a, b]
+        x = (b + a) / 2.0_dp + (b - a) * t / 2.0_dp
+
+        ! Compute function values at the transformed points
+        do i = 1, n
+            fx(i) = func(x(i))
+        end do
+
+        ! Apply the Gaussian-Legendre quadrature formula
+        integral_result = sum(w * fx) * (b - a) / 2.0_dp
+
+        ! Deallocate fx array
+        deallocate(fx)
+
+    end subroutine gauss_legendre_quadrature
+
+    ! Subroutine to initialize Gauss-Legendre nodes and weights
+    subroutine gauss_legendre_weights(t, w, n)
+        implicit none
+        integer, intent(in) :: n
+        real(dp), intent(out), dimension(n) :: t, w  !! Nodes (t) and weights (w)
+
+        ! Predefined nodes and weights for different values of n
+        select case(n)
+        case (1)
+            t = [0.0_dp]        !! Single node at the center for n = 1
+            w = [2.0_dp]        !! Weight of 2 for the single point
+        case (2)
+            t = [-0.5773502692_dp, 0.5773502692_dp]     !! Symmetric nodes for n = 2
+            w = [1.0_dp, 1.0_dp]                        !! Equal weights for n = 2
+        case (3)
+            t = [-0.7745966692_dp, 0.0_dp, 0.7745966692_dp]             !! Symmetric nodes for n = 3
+            w = [0.5555555556_dp, 0.8888888889_dp, 0.5555555556_dp]     !! Weights for n = 3
+        case (4)
+            t = [-0.8611363116_dp, -0.3399810436_dp, 0.3399810436_dp, 0.8611363116_dp]      !! Nodes for n = 4
+            w = [0.3478548451_dp, 0.6521451549_dp, 0.6521451549_dp, 0.3478548451_dp]        !! Weights for n = 4
+        case (5)
+            t = [-0.9061798459_dp, -0.5384693101_dp, 0.0_dp, 0.5384693101_dp, 0.9061798459_dp]          !! Nodes for n = 5
+            w = [0.2369268851_dp, 0.4786286705_dp, 0.5688888889_dp, 0.4786286705_dp, 0.2369268851_dp]   !! Weights for n = 5
+            ! You can add more cases to support higher values of n.
+        case default
+            print *, 'Gauss-Legendre quadrature for n > 5 is not implemented.'
+        end select
+
+    end subroutine gauss_legendre_weights
+
+end module gaussian_legendre_quadrature

modules/maths/numerical_integration/midpoint.f90

@@ -0,0 +1,65 @@
+!> Midpoint rule Module
+!!
+!! This module implements Midpoint rule for numerical integration.
+!!
+!! The midpoint rule approximates the integral by calculating the function
+!! value at the midpoint of each subinterval and summing these values, multiplied
+!! by the width of the subintervals.
+!!
+!! Note: This implementation is valid for one-dimensional functions
+!!
+!! Input:
+!! - `a`: Lower bound of integration (real(dp))
+!! - `b`: Upper bound of integration (real(dp))
+!! - `n`: Number of panels (integer)
+!! - `func`: The function to integrate (interface)
+!!
+!! Output:
+!! - `integral_result`: Approximate value of the integral (real(dp))
+
+module midpoint_rule
+    implicit none
+    integer, parameter :: dp = kind(0.d0)       !! Double precision parameter
+
+contains
+
+    subroutine midpoint(integral_result, a, b, n, func)
+        implicit none
+        integer, intent(in) :: n
+        real(dp), intent(in) :: a, b
+        real(dp), intent(out) :: integral_result
+
+        real(dp), dimension(:), allocatable :: x, fx
+        real(dp) :: h
+        integer :: i
+
+        ! Interface for the function
+        interface
+            real(kind(0.d0)) function func(x)
+                real(kind(0.d0)), intent(in) :: x
+            end function func
+        end interface
+
+        ! Step size
+        h = (b - a) / (1.0_dp*n)
+
+        ! Allocate array for midpoints
+        allocate(x(1:n), fx(1:n))
+
+        ! Calculate midpoints
+        x = [(a + (i - 0.5_dp) * h, i = 1, n)]
+
+        ! Apply function to each midpoint
+        do i = 1, n
+            fx(i) = func(x(i))
+        end do
+
+        ! Final integral value
+        integral_result = h * sum(fx)
+
+        ! Deallocate arrays
+        deallocate(x, fx)
+
+    end subroutine midpoint
+
+end module midpoint_rule