Fortran icon indicating copy to clipboard operation
Fortran copied to clipboard

Implemented Several Numerical Integration Algorithms in `maths/`

Open Ramy-Badr-Ahmed opened this issue 1 year ago • 2 comments

Description:

This pull request introduces the implementation of several numerical integration methods for approximating definite integrals over a specified interval.

Definition:

The following numerical integration methods are included:

  • Trapezoidal Rule (trapezoid.f90): Approximates the integral by dividing the interval into trapezoids and summing their areas.

  • Simpson's Rule (simpson.f90): Approximates the integral by fitting parabolas to segments of the function and is typically more accurate than the Trapezoidal Rule

  • Midpoint Rule (midpoint.f90): Approximates the integral by evaluating the function at the midpoint of each subinterval.

  • Monte Carlo Integration (with OpenMP parallelization, monte_carlo.f90): A probabilistic method that estimates the integral by sampling random points in the interval and averaging the function values. It utilizes the omp_lib with the !$omp parallel do directive for parallel execution.

  • Gaussian-Legendre Quadrature (using precomputed table values, gauss_legendre.f90): A highly accurate method that approximates the integral using a weighted sum of function values at specified points (nodes) within the integration interval. Implements Gauss-Legendre quadrature using precomputed nodes and weights for orders n=1 through n=5.

Implementation Details

  • An interface has been implemented in each module method to accept any user-defined mathematical function for integration, where functions can be passed as arguments.
  • Each implementation follows the structured module approach, including interfaces for passing the integration limits, number of points/samples, and the function to integrate.
  • The implementations leverage array slicing and vectorized operations to reduce manual loops.

Example Usage

Examples for each method are provided, demonstrating how to call the respective subroutines to approximate a definite integral of this sample function:

Input:
f(x) = exp(-x^2) * cos(2.0_dp * x)
a = -1.0_dp    !! Lower bound
b = 1.0_dp     !! Upper bound

examples/../trapezoidal.f90:

Input:
n = 1E6    !! Number of subdivisions

Output: Trapezoidal rule yields:     0.858195

examples/../simpson.f90:

Input:
n = 100    !! Number of subdivisions (must be even)

Output: Simpson's rule yields:   0.85819555

examples/../midpoint.f90:

Input:
n = 400    !! Number of subdivisionss

Output: Midpoint rule yields:     0.858196

examples/../monte_carlo.f90:

Input:
n = 1E6   !! Number of random samples

Output: Monte Carlo result:     0.857180 +-     0.000794

examples/../gaussian_legendre.f90:

Input:
n = 5    !! Number of quadrature points

Output: Gaussian Quadrature result:     0.858574

Ramy-Badr-Ahmed avatar Sep 24 '24 11:09 Ramy-Badr-Ahmed

@SatinWukerORIG

All CI checks have passed 👍

Ramy-Badr-Ahmed avatar Sep 27 '24 07:09 Ramy-Badr-Ahmed

It looks great. A lot of changes lol. Although the contribution guideline suggests not using variable names like a, b, and i, it is understandable in this case because it is for math modules. Do you think it is more readable if you change these variable names to full name?

It depends 🙂

I have renamed them in the examples files, but left them shortened in the modules files as there are explainations in the input section of the doc blocs for what they stand for. I think the forumlae in the module files would better stay in compact form with single-letter variables where comments will be helpful in this case.

Ramy-Badr-Ahmed avatar Sep 30 '24 07:09 Ramy-Badr-Ahmed

Hi @SatinWukerORIG Anything here left for this PR?

Ramy-Badr-Ahmed avatar Oct 06 '24 17:10 Ramy-Badr-Ahmed

Note

PR is complemented with #32

Ramy-Badr-Ahmed avatar Oct 12 '24 12:10 Ramy-Badr-Ahmed