Merge branch 'feature/numerical_integration'

Ramy-Badr-Ahmed · Sep 30, 2024 · 7d6bbe7 · 7d6bbe7
2 parents e343f24 + 5645569
commit 7d6bbe7
Show file tree

Hide file tree

Showing 5 changed files with 43 additions and 42 deletions.
diff --git a/examples/maths/numerical_integration/gaussian_legendre.f90 b/examples/maths/numerical_integration/gaussian_legendre.f90
@@ -19,27 +19,27 @@ program example_gaussian_quadrature
     use gaussian_legendre_quadrature
     implicit none
 
-    real(dp) :: a, b, integral_result
-    integer :: n
+    real(dp) :: lower_bound, upper_bound, integral_result
+    integer :: quadrature_points_number
 
     ! Set the integration limits and number of quadrature points
-    a = -1.0_dp
-    b = 1.0_dp
-    n = 5       !! Number of quadrature points
+    lower_bound = -1.0_dp
+    upper_bound = 1.0_dp
+    quadrature_points_number = 5       !! Number of quadrature points (order of accuracy) up to 5
 
-    ! Call Gaussian quadrature to compute the integral
-    call gauss_legendre_quadrature(integral_result, a, b, n, func)
+    ! Call Gaussian quadrature to compute the integral with the function passed as an argument
+    call gauss_legendre_quadrature(integral_result, lower_bound, upper_bound, quadrature_points_number, function)
 
     write (*, '(A, F12.6)') "Gaussian Quadrature result: ", integral_result          !! ≈ 0.858574
 
 contains
 
-    function func(x) result(fx)
+    function function(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 function function
 
 end program example_gaussian_quadrature
diff --git a/examples/maths/numerical_integration/midpoint.f90 b/examples/maths/numerical_integration/midpoint.f90
@@ -18,27 +18,28 @@
 program example_midpoint
     use midpoint_rule
     implicit none
-    real(dp) :: a, b, integral_result
-    integer :: n
+
+    real(dp) :: lower_bound, upper_bound, integral_result
+    integer :: panels_number
 
     ! Set the integration limits and number of panels
-    a = -1.0_dp
-    b = 1.0_dp
-    n = 400     !! Number of subdivisions
+    lower_bound = -1.0_dp
+    upper_bound = 1.0_dp
+    panels_number = 400     !! Number of subdivisions
 
     ! Call the midpoint rule subroutine with the function passed as an argument
-    call midpoint(integral_result, a, b, n, func)
+    call midpoint(integral_result, lower_bound, upper_bound, panels_number, function)
 
     write (*, '(A, F12.6)') "Midpoint rule yields: ", integral_result  !! ≈ 0.858196
 
 contains
 
-    function func(x) result(fx)
+    function function(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 function function
 
 end program example_midpoint
diff --git a/examples/maths/numerical_integration/monte_carlo.f90 b/examples/maths/numerical_integration/monte_carlo.f90
@@ -19,27 +19,27 @@ program example_monte_carlo
     use monte_carlo_integration
     implicit none
 
-    real(dp) :: a, b, integral_result, error_estimate
-    integer :: n
+    real(dp) :: lower_bound, upper_bound, integral_result, error_estimate
+    integer :: random_samples_number
 
     ! Set the integration limits and number of random samples
-    a = -1.0_dp
-    b = 1.0_dp
-    n = 1000000     !! 1E6 Number of random samples
+    lower_bound = -1.0_dp
+    upper_bound = 1.0_dp
+    random_samples_number = 1000000     !! 1E6 Number of random samples
 
-    ! Call Monte Carlo integration
-    call monte_carlo(integral_result, error_estimate, a, b, n, func)
+    ! Call Monte Carlo integration with the function passed as an argument
+    call monte_carlo(integral_result, error_estimate, lower_bound, upper_bound, random_samples_number, function)
 
     write (*, '(A, F12.6, A, F12.6)') "Monte Carlo result: ", integral_result, " +- ", error_estimate     !! ≈ 0.858421
 
 contains
 
-    function func(x) result(fx)
+    function function(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 function function
 
 end program example_monte_carlo
diff --git a/examples/maths/numerical_integration/simpson.f90 b/examples/maths/numerical_integration/simpson.f90
@@ -19,27 +19,27 @@ program example_simpson
     use simpson_rule
     implicit none
 
-    real(dp) :: a, b, integral_result
-    integer :: n
+    real(dp) :: lower_bound, upper_bound, integral_result
+    integer :: panels_number
 
     ! Set the integration limits and number of panels
-    a = -1.0_dp
-    b = 1.0_dp
-    n = 100     !! Number of subdivisions (must be even)
+    lower_bound = -1.0_dp
+    upper_bound = 1.0_dp
+    panels_number = 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)
+    call simpson(integral_result, lower_bound, upper_bound, panels_number, function)
 
     write (*, '(A, F12.8)') "Simpson's rule yields: ", integral_result       !! ≈ 0.85819555
 
 contains
 
-    function func(x) result(fx)
+    function function(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 function function
 
 end program example_simpson
diff --git a/examples/maths/numerical_integration/trapezoid.f90 b/examples/maths/numerical_integration/trapezoid.f90
@@ -19,27 +19,27 @@ program example_tapezoid
     use trapezoidal_rule
     implicit none
 
-    real(dp) :: a, b, integral_result
-    integer :: n
+    real(dp) :: lower_bound, upper_bound, integral_result
+    integer :: panels_number
 
     ! Set the integration limits and number of panels
-    a = -1.0_dp
-    b = 1.0_dp
-    n = 1000000     !! 1E6 Number of subdivisions
+    lower_bound = -1.0_dp
+    upper_bound = 1.0_dp
+    panels_number = 1000000     !! 1E6 Number of subdivisions
 
     ! Call the trapezoidal rule with the function passed as an argument
-    call trapezoid(integral_result, a, b, n, func)
+    call trapezoid(integral_result, lower_bound, upper_bound, panels_number, function)
 
     write (*, '(A, F12.6)') 'Trapezoidal rule yields: ', integral_result     !! ≈ 0.858195
 
 contains
 
-    function func(x) result(fx)
+    function function(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 function function
 
 end program example_tapezoid