How to do multi-dimension integration with BigFloat?

New to Julia
1

Hello, I want to calculate a multiple dimension integration,
and I have to use high precision float.

I have tried HCubature, which is a pure-julia implementation of integration.
But when I use BigFloat as input, it will report an error:

ERROR: LoadError: setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.

What does that mean?
Is there any other packages that I can use?

2

It is recommended that you show a script that reproduces the problem.
It appears to work with BigFloat for simple cases.

julia> using HCubature

julia> f(x, y) = x^2 + 5y^2
f (generic function with 1 method)

julia> f(v) = f(v...)
f (generic function with 2 methods)

julia> a0, b0 = big"0.0", big"1.0"
(0.0, 1.0)

julia> a1, b1 = big"0.0", big"2.0"
(0.0, 2.0)

julia> hcubature(f, (a0, a1), (b0, b1))
(14.0, 0.0)

Also, do you need arbitrary precision?
If you just need more precision than Float64, Float128 from Quadmath.jl might be useful.

3

Although HCubature.jl does indeed support arbitrary precision (its multidimensional quadrature rule is carefully written to use the precision of the endpoints), it corresponds to a low-order quadrature rule that will be very slow (require a lot of subdivisions) to reach high accuracies.

If you are using BigFloat or similar types to obtain high accuracy (\gg 16 digits), then you really want to use a high-order quadrature rule (computed in BigFloat precision too). The QuadGK.jl package supports this, as discussed in this tutorial example and this post, for example. You can use nested quadgk calls to do multidimensional integration, and there is a package IteratedIntegration.jl wrapping QuadGK that does this more efficiently.

However, if you want high accuracy arbitrary precision integration, you should also have an integrand that is very smooth (or for which any singularities are built into the quadrature rule). In that case, you might not want an adaptive routine at all, but simply use a tensor product of 1d quadrature rules (computed with QuadGK to arbitrary precision), e.g. see the example in this post.

2 Likes