trapz - integrate sampled values using trapezoidal ruleExperimental
Returns the trapezoidal rule integral of an array y representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing dx or arbitrary abscissas x.
result = trapz (y, x)
result = trapz (y, dx)
y: Shall be a rank-one array of type real.
x: Shall be a rank-one array of type real having the same kind and size as y.
dx: Shall be a scalar of type real having the same kind as y.
The result is a scalar of type real having the same kind as y.
If the size of y is zero or one, the result is zero.
program example_trapz
use stdlib_quadrature, only: trapz
implicit none
real, parameter :: x(5) = [0., 1., 2., 3., 4.]
real :: y(5) = x**2
print *, trapz(y, x)
! 22.0
print *, trapz(y, 0.5)
! 11.0
end program example_trapz
trapz_weights - trapezoidal rule weights for given abscissasExperimental
Given an array of abscissas x, computes the array of weights w such that if y represented function values tabulated at x, then sum(w*y) produces a trapezoidal rule approximation to the integral.
result = trapz_weights (x)
x: Shall be a rank-one array of type real.
The result is a real array with the same size and kind as x.
If the size of x is one, then the sole element of the result is zero.
program example_trapz_weights
use stdlib_quadrature, only: trapz_weights
implicit none
real, parameter :: x(5) = [0., 1., 2., 3., 4.]
real :: y(5) = x**2
real :: w(5)
w = trapz_weights(x)
print *, sum(w*y)
! 22.0
end program example_trapz_weights
simps - integrate sampled values using Simpson's ruleExperimental
Returns the Simpson's rule integral of an array y representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing dx or arbitrary abscissas x.
Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of even. If even is negative (positive), the 3/8 rule is used at the beginning (end) of the array. If even is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and these two results were averaged.
result = simps (y, x [, even])
result = simps (y, dx [, even])
y: Shall be a rank-one array of type real.
x: Shall be a rank-one array of type real having the same kind and size as y.
dx: Shall be a scalar of type real having the same kind as y.
even: (Optional) Shall be a default-kind integer.
The result is a scalar of type real having the same kind as y.
If the size of y is zero or one, the result is zero.
If the size of y is two, the result is the same as if trapz had been called instead.
program example_simps
use stdlib_quadrature, only: simps
implicit none
real, parameter :: x(5) = [0., 1., 2., 3., 4.]
real :: y(5) = 3.*x**2
print *, simps(y, x)
! 64.0
print *, simps(y, 0.5)
! 32.0
end program example_simps
simps_weights - Simpson's rule weights for given abscissasExperimental
Given an array of abscissas x, computes the array of weights w such that if y represented function values tabulated at x, then sum(w*y) produces a Simpson's rule approximation to the integral.
Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of even. If even is negative (positive), the 3/8 rule is used at the beginning (end) of the array and the 1/3 rule used elsewhere. If even is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and then these two results were averaged.
result = simps_weights (x [, even])
x: Shall be a rank-one array of type real.
even: (Optional) Shall be a default-kind integer.
The result is a real array with the same size and kind as x.
If the size of x is one, then the sole element of the result is zero.
If the size of x is two, then the result is the same as if trapz_weights had been called instead.
program example_simps_weights
use stdlib_quadrature, only: simps_weights
implicit none
real, parameter :: x(5) = [0., 1., 2., 3., 4.]
real :: y(5) = 3.*x**2
real :: w(5)
w = simps_weights(x)
print *, sum(w*y)
! 64.0
end program example_simps_weights
gauss_legendre - Gauss-Legendre quadrature (a.k.a. Gaussian quadrature) nodes and weightsExperimental
Computes Gauss-Legendre quadrature (also known as simply Gaussian quadrature) nodes and weights,
for any N (number of nodes).
Using the nodes x and weights w, you can compute the integral of some function f as follows:
integral = sum(f(x) * w).
Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself. Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision (maximum difference from those values is 2 epsilon).
subroutine gauss_legendre (x, w[, interval])
x: Shall be a rank-one array of type real(real64). It is an output argument, representing the quadrature nodes.
w: Shall be a rank-one array of type real(real64), with the same dimension as x.
It is an output argument, representing the quadrature weights.
interval: (Optional) Shall be a two-element array of type real(real64).
If present, the nodes and weigts are calculated for integration from interval(1) to interval(2).
If not specified, the default integral is -1 to 1.
program example_gauss_legendre
use iso_fortran_env, dp => real64
use stdlib_quadrature, only: gauss_legendre
implicit none
integer, parameter :: N = 6
real(dp), dimension(N) :: x, w
call gauss_legendre(x, w)
print *, "integral of x**2 from -1 to 1 is", sum(x**2*w)
end program example_gauss_legendre
gauss_legendre_lobatto - Gauss-Legendre-Lobatto quadrature nodes and weightsExperimental
Computes Gauss-Legendre-Lobatto quadrature nodes and weights,
for any N (number of nodes).
Using the nodes x and weights w, you can compute the integral of some function f as follows:
integral = sum(f(x) * w).
Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself. Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision (maximum difference from those values is 2 epsilon).
subroutine gauss_legendre_lobatto (x, w[, interval])
x: Shall be a rank-one array of type real(real64). It is an output argument, representing the quadrature nodes.
w: Shall be a rank-one array of type real(real64), with the same dimension as x.
It is an output argument, representing the quadrature weights.
interval: (Optional) Shall be a two-element array of type real(real64).
If present, the nodes and weigts are calculated for integration from interval(1) to interval(2).
If not specified, the default integral is -1 to 1.
program example_gauss_legendre_lobatto
use iso_fortran_env, dp => real64
use stdlib_quadrature, only: gauss_legendre_lobatto
implicit none
integer, parameter :: N = 6
real(dp), dimension(N) :: x, w
call gauss_legendre_lobatto(x, w)
print *, "integral of x**2 from -1 to 1 is", sum(x**2*w)
end program example_gauss_legendre_lobatto