Skip to content

Commit

Permalink
Merge branch 'feature/numerical_integration'
Browse files Browse the repository at this point in the history
  • Loading branch information
Ramy-Badr-Ahmed committed Sep 30, 2024
2 parents e343f24 + 5645569 commit 7d6bbe7
Show file tree
Hide file tree
Showing 5 changed files with 43 additions and 42 deletions.
18 changes: 9 additions & 9 deletions examples/maths/numerical_integration/gaussian_legendre.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
17 changes: 9 additions & 8 deletions examples/maths/numerical_integration/midpoint.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
18 changes: 9 additions & 9 deletions examples/maths/numerical_integration/monte_carlo.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
16 changes: 8 additions & 8 deletions examples/maths/numerical_integration/simpson.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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