-
-
Notifications
You must be signed in to change notification settings - Fork 23
/
Copy pathtrapezoid.f90
72 lines (60 loc) · 2.1 KB
/
trapezoid.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
!> Trapezoidal Rule Module
!!
!! This module implements the Trapezoidal rule for numerical integration.
!!
!! Created by: Ramy-Badr-Ahmed (https://github.com/Ramy-Badr-Ahmed)
!! in Pull Request: #25
!! https://github.com/TheAlgorithms/Fortran/pull/25
!!
!! Please mention me (@Ramy-Badr-Ahmed) in any issue or pull request
!! addressing bugs/corrections to this file. Thank you!
!!
!! The Trapezoidal rule approximates the definite integral of a function by
!! dividing the area under the curve into trapezoids and summing their areas.
!!
!! Note: This implementation is valid for one-dimensional functions
!!
!! Input:
!! - `a`: Lower bound of the integration (real(dp))
!! - `b`: Upper bound of the integration (real(dp))
!! - `n`: Number of panels (integer)
!! - `func`: The function to integrate (interface)
!!
!! Output:
!! - `integral_result`: Approximate value of the definite integral (real(dp))
!!
module trapezoidal_rule
implicit none
integer, parameter :: dp = kind(0.d0) !! Double precision parameter
contains
! Trapezoidal rule with function passed via interface
subroutine trapezoid(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)/real(n, dp)
! Allocate arrays
allocate (x(0:n), fx(0:n))
! Create an array of x values
x = [(a + (real(i, dp))*h, i=0, n)]
! Apply the function to each x value
do i = 0, n
fx(i) = func(x(i))
end do
! Apply trapezoidal rule using array slicing
integral_result = ((fx(0) + fx(n))*0.5_dp + sum(fx(1:n)))*h
! Deallocate arrays
deallocate (x, fx)
end subroutine trapezoid
end module trapezoidal_rule