• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    迪恩网络公众号

JuliaMath/HCubature.jl: pure-Julia multidimensional h-adaptive integration

原作者: [db:作者] 来自: 网络 收藏 邀请

开源软件名称:

JuliaMath/HCubature.jl

开源软件地址:

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

开源编程语言:

Julia 100.0%

开源软件介绍:

HCubature

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

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.

Usage

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.)

Algorithm

The algorithm of hcubature is based on the one described in:

Author and Copyright

HCubature was written by Steven G. Johnson (SGJ), and is free/open-source software under the MIT/expat license.

SGJ also wrote an earlier C implementation of a similar algorithm that is also callable from Julia via the Cubature.jl package. The HCubature package is a from-scratch re-implementation, not a translation, of this code, both to take advantage of unique features of Julia and to eliminate licensing restrictions arising from the use of C code taken from the HIntLib library. (In both cases, the original DCUHRE Fortran code of Genz was not examined, only the mathematical description in the papers.)




鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap