1661237220

This package provides support for one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature. The code was originally part of Base Julia. It supports integration of arbitrary numeric types, including arbitrary precision (`BigFloat`

), and even integration of arbitrary normed vector spaces (e.g. matrix-valued integrands).

The package provides three functions: `quadgk`

, `gauss`

, and `kronrod`

. `quadgk`

performs the integration, `gauss`

computes Gaussian quadrature points and weights for integrating over the interval [a, b], and `kronrod`

computes Kronrod points, weights, and embedded Gaussian quadrature weights for integrating over [-1, 1]. Typical usage looks like:

```
using QuadGK
integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
```

which computes the integral of exp(–x²) from x=0 to x=1 to a relative tolerance of 10⁻⁸, and returns the approximate `integral = 0.746824132812427`

and error estimate `err = 7.887024366937112e-13`

(which is actually smaller than the requested tolerance: convergence was very rapid because the integrand is smooth).

For more information, see the documentation.

For integrands whose values are *small* arrays whose length is known at compile-time, it is usually most efficient to modify your integrand to return an `SVector`

from the StaticArrays.jl package.

However, for integrands that return large or variabley-length arrays, we also provide a function `quadgk!(f!, result, a,b...)`

in order to exploit in-place operations where possible. The `result`

argument is used to store the estimated integral `I`

in-place, and the integrand function is now of the form `f!(r, x)`

and should write `f(x)`

in-place into the result array `r`

.

If you are computing many similar integrals of smooth functions, you may not need an adaptive integration — with a little experimentation, you may be able to decide on an appropriate number `N`

of integration points in advance, and re-use this for all of your integrals. In this case you can use `x, w = gauss(N, a, b)`

to find the quadrature points `x`

and weights `w`

, so that `sum(f.(x) .* w)`

is an `N`

-point approximation to `∫f(x)dx`

from `a`

to `b`

.

For computing many integrands of similar functions with *singularities*, `x, w = gauss(W, N, a, b)`

function allows you to pass a *weight function* `W(x)`

as the first argument, so that `sum(f.(x) .* w)`

is an `N`

-point approximation to `∫W(x)f(x)dx`

from `a`

to `b`

. In this way, you can put all of the singularities etcetera into `W`

and precompute an accurate quadrature rule as long as the remaining `f(x)`

terms are smooth. For example,

```
using QuadGK
x, w = gauss(x -> exp(-x) / sqrt(x), 10, 0, -log(1e-10), rtol=1e-9)
```

computes the points and weights for performing `∫exp(-x)f(x)/√x dx`

integrals from `0`

to `-log(1e-10) ≈ 23`

, so that there is a `1/√x`

singularity in the integrand at `x=0`

and a rapid decay for increasing `x`

. (The `gauss`

function currently does not support infinite integration intervals, but for a rapidly decaying weight function you can approximate an infinite interval to any desired accuracy by a sufficiently broad interval, with a tradeoff in computational expense.) For example, with `f(x) = sin(x)`

, the exact answer is `0.570370556005742…`

. Using the points and weights above with `sum(sin.(x) .* w)`

, we obtain `0.5703706212868831`

, which is correct to 6–7 digits using only 10 `f(x)`

evaluations. Obtaining similar accuracy for the same integral from `quadgk`

requires nearly 300 function evaluations. However, the `gauss`

function itself computes many (`2N`

) numerical integrals of your weight function (multiplied by polynomials), so this is only more efficient if your `f(x)`

is very expensive or if you need to compute a large number of integrals with the same `W`

.

See the `gauss`

documentation for more information. See also our example using a weight function interpolated from tabulated data.

The FastGaussQuadrature.jl package provides non-adaptive Gaussian quadrature variety of built-in weight functions — it is a good choice you need to go to very high orders N, e.g. to integrate rapidly oscillating functions, or use weight functions that incorporate some standard singularity in your integrand. QuadGK, on the other hand, keeps the order N of the quadrature rule fixed and improves accuracy by subdividing the integration domain, which can be better if fine resolution is required only in a part of your domain (e.g if your integrand has a sharp peak or singularity somewhere that is not known in advance).

For multidimensional integration, see the HCubature.jl, Cubature.jl, and Cuba.jl packages.

Documentation:

Author: JuliaMath

Source Code: https://github.com/JuliaMath/QuadGK.jl

License: MIT license

1661237220

This package provides support for one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature. The code was originally part of Base Julia. It supports integration of arbitrary numeric types, including arbitrary precision (`BigFloat`

), and even integration of arbitrary normed vector spaces (e.g. matrix-valued integrands).

The package provides three functions: `quadgk`

, `gauss`

, and `kronrod`

. `quadgk`

performs the integration, `gauss`

computes Gaussian quadrature points and weights for integrating over the interval [a, b], and `kronrod`

computes Kronrod points, weights, and embedded Gaussian quadrature weights for integrating over [-1, 1]. Typical usage looks like:

```
using QuadGK
integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
```

which computes the integral of exp(–x²) from x=0 to x=1 to a relative tolerance of 10⁻⁸, and returns the approximate `integral = 0.746824132812427`

and error estimate `err = 7.887024366937112e-13`

(which is actually smaller than the requested tolerance: convergence was very rapid because the integrand is smooth).

For more information, see the documentation.

For integrands whose values are *small* arrays whose length is known at compile-time, it is usually most efficient to modify your integrand to return an `SVector`

from the StaticArrays.jl package.

However, for integrands that return large or variabley-length arrays, we also provide a function `quadgk!(f!, result, a,b...)`

in order to exploit in-place operations where possible. The `result`

argument is used to store the estimated integral `I`

in-place, and the integrand function is now of the form `f!(r, x)`

and should write `f(x)`

in-place into the result array `r`

.

If you are computing many similar integrals of smooth functions, you may not need an adaptive integration — with a little experimentation, you may be able to decide on an appropriate number `N`

of integration points in advance, and re-use this for all of your integrals. In this case you can use `x, w = gauss(N, a, b)`

to find the quadrature points `x`

and weights `w`

, so that `sum(f.(x) .* w)`

is an `N`

-point approximation to `∫f(x)dx`

from `a`

to `b`

.

For computing many integrands of similar functions with *singularities*, `x, w = gauss(W, N, a, b)`

function allows you to pass a *weight function* `W(x)`

as the first argument, so that `sum(f.(x) .* w)`

is an `N`

-point approximation to `∫W(x)f(x)dx`

from `a`

to `b`

. In this way, you can put all of the singularities etcetera into `W`

and precompute an accurate quadrature rule as long as the remaining `f(x)`

terms are smooth. For example,

```
using QuadGK
x, w = gauss(x -> exp(-x) / sqrt(x), 10, 0, -log(1e-10), rtol=1e-9)
```

computes the points and weights for performing `∫exp(-x)f(x)/√x dx`

integrals from `0`

to `-log(1e-10) ≈ 23`

, so that there is a `1/√x`

singularity in the integrand at `x=0`

and a rapid decay for increasing `x`

. (The `gauss`

function currently does not support infinite integration intervals, but for a rapidly decaying weight function you can approximate an infinite interval to any desired accuracy by a sufficiently broad interval, with a tradeoff in computational expense.) For example, with `f(x) = sin(x)`

, the exact answer is `0.570370556005742…`

. Using the points and weights above with `sum(sin.(x) .* w)`

, we obtain `0.5703706212868831`

, which is correct to 6–7 digits using only 10 `f(x)`

evaluations. Obtaining similar accuracy for the same integral from `quadgk`

requires nearly 300 function evaluations. However, the `gauss`

function itself computes many (`2N`

) numerical integrals of your weight function (multiplied by polynomials), so this is only more efficient if your `f(x)`

is very expensive or if you need to compute a large number of integrals with the same `W`

.

See the `gauss`

documentation for more information. See also our example using a weight function interpolated from tabulated data.

The FastGaussQuadrature.jl package provides non-adaptive Gaussian quadrature variety of built-in weight functions — it is a good choice you need to go to very high orders N, e.g. to integrate rapidly oscillating functions, or use weight functions that incorporate some standard singularity in your integrand. QuadGK, on the other hand, keeps the order N of the quadrature rule fixed and improves accuracy by subdividing the integration domain, which can be better if fine resolution is required only in a part of your domain (e.g if your integrand has a sharp peak or singularity somewhere that is not known in advance).

For multidimensional integration, see the HCubature.jl, Cubature.jl, and Cuba.jl packages.

Documentation:

Author: JuliaMath

Source Code: https://github.com/JuliaMath/QuadGK.jl

License: MIT license

1661256480

This is library intended to provided multidimensional numerical integration routines in pure Julia

For the time being this library can only perform integrals in three dimensions.

**TODO**

- Add rules for other dimensions
- Make sure it works properly with complex valued functions
- Parallelize
- Improve the error estimates (the Cuba library and consequently Cuba.jl seem to calculate tighter errors)

`NIntegration.jl`

should work on Julia 1.0 and later versions and can be installed from a Julia session by running

```
julia> using Pkg
julia> Pkg.add(PackageSpec(url = "https://github.com/pabloferz/NIntegration.jl.git"))
```

Once installed, run

```
using NIntegration
```

To integrate a function `f(x, y, z)`

on the hyperrectangle defined by `xmin`

and `xmax`

, just call

```
nintegrate(
f::Function, xmin::NTuple{N}, xmax::NTuple{N};
reltol = 1e-6, abstol = eps(), maxevals = 1000000
)
```

The above returns a tuple `(I, E, n, R)`

of the calculated integral `I`

, the estimated error `E`

, the number of integrand evaluations `n`

, and a list `R`

of the subregions in which the integration domain was subdivided.

If you need to evaluate multiple functions `(f₁, f₂, ...)`

on the same integration domain, you can evaluate the function `f`

with more "features" and use its subregions list to estimate the integral for the rest of the functions in the list, e.g.

```
(I, E, n, R) = nintegrate(f, xmin, xmin)
I₁ = nintegrate(f₁, R)
```

The integration algorithm is based on the one decribed in:

- J. Berntsen, T. O. Espelid, and A. Genz, "An Adaptive Algorithm for the Approximate Calculation of Multiple Integrals,"
*ACM Trans. Math. Soft.*, 17 (4), 437-451 (1991).

The author expresses his gratitude to Professor Alan Genz for some useful pointers.

This work was financially supported by CONACYT through grant 354884.

Author: Pabloferz

Source Code: https://github.com/pabloferz/NIntegration.jl

License: MIT license

1661244900

The HCubature module is a pure-Julia implementation of multidimensional "h-adaptive" integration. That is, given an n-dimensional integral

then `hcubature(f, a, b)`

computes the integral, adaptively subdividing the integration volume into smaller and smaller pieces until convergence is achieved to the desired tolerance (specified by optional `rtol`

and `atol`

keyword arguments, described in more detail below.

Because `hcubature`

is written purely in Julia, the integrand `f(x)`

can return any vector-like object (technically, any type supporting `+`

, `-`

, `*`

real, and `norm`

: a Banach space). You can integrate real, complex, and matrix-valued integrands, for example.

Assuming you've installed the HCubature package (via `Pkg.add`

) and loaded it with `using HCubature`

, you can then use it by calling the `hcubature`

function:

`hcubature`

```
hcubature(f, a, b; norm=norm, rtol=sqrt(eps), atol=0, maxevals=typemax(Int), initdiv=1)
```

This computes the n-dimensional integral of f(x), where `n == length(a) == length(b)`

, over the hypercube whose corners are given by the vectors (or tuples) `a`

and `b`

. That is, dimension `x[i]`

is integrated from `a[i]`

to `b[i]`

. The return value of `hcubature`

is a tuple `(I, E)`

of the estimated integral `I`

and an estimated error `E`

.

`f`

should be a function `f(x)`

that takes an n-dimensional vector `x`

and returns the integrand at `x`

. The integrand can be any type that supports `+`

, `-`

, `*`

real, and `norm`

functions. For example, the integrand can be real or complex numbers, vectors, matrices, etcetera. (For performance, the StaticArrays package is recommended for use with vector/matrix-valued integrands.)

The integrand `f(x)`

will be always be passed an `SVector{n,T}`

, where `SVector`

is an efficient vector type defined in the StaticArrays package and `T`

is a floating-point type determined by promoting the endpoint `a`

and `b`

coordinates to a floating-point type. (Your integrand `f`

should be type-stable: it should always return a value of the same type, given this type of `x`

.)

The integrand will never be evaluated exactly at the boundaries of the integration volume. (So, for example, it is possible to have an integrand that blows up at the boundaries, as long as the integral is finite, though such singularities will slow convergence.)

The integration volume is adaptively subdivided, using a cubature rule due to Genz and Malik (1980), until the estimated error `E`

satisfies `E ≤ max(rtol*norm(I), atol)`

, i.e. `rtol`

and `atol`

are the relative and absolute tolerances requested, respectively. It also stops if the number of `f`

evaluations exceeds `maxevals`

. If neither `atol`

nor `rtol`

are specified, the default `rtol`

is the square root of the precision `eps(T)`

of the coordinate type `T`

described above. Initially, the volume is divided into `initdiv`

segments along each dimension.

The error is estimated by `norm(I - I′)`

, where `I′`

is an alternative estimated integral (via an "embedded" lower-order cubature rule.) By default, the `norm`

function used (for both this and the convergence test above) is `norm`

, but you can pass an alternative norm by the `norm`

keyword argument. (This is especially useful when `f`

returns a vector of integrands with different scalings.)

`hquadrature`

```
hquadrature(f, a, b; norm=norm, rtol=sqrt(eps), atol=0, maxevals=typemax(Int), initdiv=1)
```

Compute the (1d) integral of f(x) from `a`

to `b`

. The return value of `hcubature`

is a tuple `(I, E)`

of the estimated integral `I`

and an estimated error `E`

.

The other parameters are the same as `hcubature`

(above). `hquadrature`

is just a convenience wrapper around `hcubature`

so that you can work with scalar `x`

, `a`

, and `b`

, rather than 1-component vectors.

Alternatively, for 1d integrals you can import the QuadGK module and call the `quadgk`

function, which provides additional flexibility e.g. in choosing the order of the quadrature rule. (`QuadGK`

is used internally anyway by `HCubature`

to compute the quadrature rule.)

The algorithm of `hcubature`

is based on the one described in:

- A. C. Genz and A. A. Malik, "An adaptive algorithm for numeric integration over an N-dimensional rectangular region,"
*J. Comput. Appl. Math.*, vol. 6 (no. 4), 295-302 (1980).

Author: JuliaMath

Source Code: https://github.com/JuliaMath/HCubature.jl

License: View license

1665717240

AMD clBLAS bindings for Julia.

if the download of the binary fails on linux, you may also try to install clblas with `sudo apt-get install libclblas-dev`

and then rerun `Pkg.build("CLBLAS")`

.

- Low-level bindings
- Partial implementation of high-level API similar to
`Base.LinAlg.BLAS`

```
using OpenCL
import OpenCL.cl.CLArray
import CLBLAS
const clblas = CLBLAS
clblas.setup()
device, ctx, queue = clblas.get_next_compute_context()
alpha = 1.0
beta = 0.0
hA = rand(5, 10)
hB = rand(10, 5)
A = CLArray(queue, hA)
B = CLArray(queue, hB)
C = cl.zeros(queue, 5, 5)
clblas.gemm!('N', 'N', alpha, A, B, beta, C)
hC = cl.to_host(C)
if isapprox(hC, hA * hB)
info("Success!")
else
error("Results diverged")
end
```

`Complex64`

/`CL_float2`

doesn't work by default. This is caused by some weird incompatibility between clang (default for Julia) and gcc (default for clBLAS), so the only way to fix it right now is to manually compile clBLAS using clang,

Author: JuliaGPU

Source Code: https://github.com/JuliaGPU/CLBLAS.jl

License: Apache-2.0 license

1661252700

`Cuba.jl`

is a library for multidimensional numerical integration with different algorithms in Julia.

This is just a Julia wrapper around the C Cuba library, version 4.2, by **Thomas Hahn**. All the credits goes to him for the underlying functions, blame me for any problem with the Julia interface. Feel free to report bugs and make suggestions at https://github.com/giordano/Cuba.jl/issues.

All algorithms provided by Cuba library are supported in `Cuba.jl`

:

`vegas`

(type: Monte Carlo; variance reduction with importance sampling)`suave`

(type: Monte Carlo; variance reduction with globally adaptive subdivision + importance sampling)`divonne`

(type: Monte Carlo or deterministic; variance reduction with stratified sampling, aided by methods from numerical optimization)`cuhre`

(type: deterministic; variance reduction with globally adaptive subdivision)

Integration is performed on the n-dimensional unit hypercube [0, 1]^n. For more details on the algorithms see the manual included in Cuba library and available in `deps/usr/share/cuba.pdf`

after successful installation of `Cuba.jl`

.

`Cuba.jl`

is available for GNU/Linux, FreeBSD, Mac OS, and Windows (`i686`

and `x86_64`

architectures).

The latest version of `Cuba.jl`

is available for Julia 1.3 and later versions, and can be installed with Julia built-in package manager. In a Julia session, after entering the package manager mode with `]`

, run the command

```
pkg> update
pkg> add Cuba
```

Older versions are also available for Julia 0.4-1.2.

After installing the package, run

```
using Cuba
```

or put this command into your Julia script.

`Cuba.jl`

provides the following functions to integrate:

```
vegas(integrand, ndim, ncomp[; keywords...])
suave(integrand, ndim, ncomp[; keywords...])
divonne(integrand, ndim, ncomp[; keywords...])
cuhre(integrand, ndim, ncomp[; keywords...])
```

These functions wrap the 64-bit integers functions provided by the Cuba library.

The only mandatory argument is:

`function`

: the name of the function to be integrated

Optional positional arguments are:

`ndim`

: the number of dimensions of the integration domain. Defaults to 1 in`vegas`

and`suave`

, to 2 in`divonne`

and`cuhre`

. Note:`ndim`

must be at least 2 with the latest two methods.`ncomp`

: the number of components of the integrand. Defaults to 1

`ndim`

and `ncomp`

arguments must appear in this order, so you cannot omit `ndim`

but not `ncomp`

. `integrand`

should be a function `integrand(x, f)`

taking two arguments:

- the input vector
`x`

of length`ndim`

- the output vector
`f`

of length`ncomp`

, used to set the value of each component of the integrand at point`x`

Also anonymous functions are allowed as `integrand`

. For those familiar with `Cubature.jl`

package, this is the same syntax used for integrating vector-valued functions.

For example, the integral

```
∫_0^1 cos(x) dx = sin(1) = 0.8414709848078965
```

can be computed with one of the following commands

```
julia> vegas((x, f) -> f[1] = cos(x[1]))
Component:
1: 0.8414910005259609 ± 5.2708169787733e-5 (prob.: 0.028607201257039333)
Integrand evaluations: 13500
Number of subregions: 0
Note: The desired accuracy was reached
julia> suave((x, f) -> f[1] = cos(x[1]))
Component:
1: 0.8411523690658836 ± 8.357995611133613e-5 (prob.: 1.0)
Integrand evaluations: 22000
Number of subregions: 22
Note: The desired accuracy was reached
julia> divonne((x, f) -> f[1] = cos(x[1]))
Component:
1: 0.841468071955942 ± 5.3955070531551656e-5 (prob.: 0.0)
Integrand evaluations: 1686
Number of subregions: 14
Note: The desired accuracy was reached
julia> cuhre((x, f) -> f[1] = cos(x[1]))
Component:
1: 0.8414709848078966 ± 2.2204460420128823e-16 (prob.: 3.443539937576958e-5)
Integrand evaluations: 195
Number of subregions: 2
Note: The desired accuracy was reached
```

The integrating functions `vegas`

, `suave`

, `divonne`

, and `cuhre`

return an `Integral`

object whose fields are

```
integral :: Vector{Float64}
error :: Vector{Float64}
probl :: Vector{Float64}
neval :: Int64
fail :: Int32
nregions :: Int32
```

The first three fields are vectors with length `ncomp`

, the last three ones are scalars. The `Integral`

object can also be iterated over like a tuple. In particular, if you assign the output of integration functions to the variable named `result`

, you can access the value of the `i`

-th component of the integral with `result[1][i]`

or `result.integral[i]`

and the associated error with `result[2][i]`

or `result.error[i]`

. The details of other quantities can be found in Cuba manual.

All other arguments listed in Cuba documentation can be passed as optional keywords.

A more detailed manual of `Cuba.jl`

, with many complete examples, is available at https://giordano.github.io/Cuba.jl/stable/.

There are other Julia packages for multidimenensional numerical integration:

Author: Giordano

Source Code: https://github.com/giordano/Cuba.jl

License: View license