- Numerical integration
- trapz - integrate sampled values using trapezoidal rule
- trapz_weights - trapezoidal rule weights for given abscissas
- simps - integrate sampled values using Simpson's rule
- simps_weights - Simpson's rule weights for given abscissas
- gauss_legendre - Gauss-Legendre quadrature (a.k.a. Gaussian quadrature) nodes and weights
- gauss_legendre_lobatto - Gauss-Legendre-Lobatto quadrature nodes and weights

`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 arbitary 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 demo_trapz
use stdlib_quadrature, only: trapz
implicit none
real :: 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 demo_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 demo_trapz_weights
use stdlib_quadrature, only: trapz_weights
implicit none
real :: 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 demo_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 arbitary 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 demo_simps
use stdlib_quadrature, only: simps
implicit none
real :: 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 demo_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 demo_simps_weights
use stdlib_quadrature, only: simps_weights
implicit none
real :: 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 demo_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 integrate
use iso_fortran_env, dp => real64
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
```

`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 integrate
use iso_fortran_env, dp => real64
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
```