Advice on approximating Complicated integral

General Usage
1

hi everyone i’m trying to solve a system of fixed point equations that involve integrals of this type.
immagine

so far nesting quadgk calls results with high probability in a domain error, the only thing that helps is using BigFloats, but then the performance degrades to the point where it becomes unthinkable to solve the system involving this integrals,

is there a good soul with some advice for me? (besides forgetting this integral ever existed :rofl:)
thank you!

2

Are your integrals from -\infty to +\infty? You could try using a specialized quadrature rule, e.g. Gauss–Hermite quadrature would be reasonable here given the \exp(-h^2/2) weight factor. The GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature package can generate such rules for you. (This is for a fixed quadrature order, so you’d have to do some experimentation to find a good quadrature order for your problem.)

3 Likes
3

Thank you Steven, your advice turned out to be a winning strategy, i’m also started on actually solving my equations.

this is the code i use to calculate those ‘gaussian’ expectation values now

using FastGaussQuadrature

struct GaussHermite{T,order}
    X::Vector{T}
    W::Vector{T}
end

function GaussHermite(T,order)
    X,W=gausshermite( order )
    X=T.(X).*sqrt(T(2))
    W=T.(W)./sqrt(T(π))
    idx=W.>=eps(T)/2
    X=X[idx]
    W=W[idx]
    GaussHermite{T,length(X)}(X,W)
end

function (g::GaussHermite{T,order})(f) where {T,order}
    I0=g.W[1].*f(g.X[1])
    @inbounds for i in 2:order
        I0=I0.+g.W[i].*f(g.X[i])
    end
    I0
end

and this is some testing code to find out what the best order actually is

function testgh()
    F(x)=cos(x*x/30)*tanh(0.5+10x)*x*x
    gmoment=0.0003549752261842731 # Calculated in Mathematica with PrecisionGoal->11
    least=Inf
    oleast=0
    for ord in 1:2000
        gh=GaussHermite(Float64,ord)
        Δ = gh(F)-gmoment
        mΔ=maximum(abs,Δ)
        mΔ<least && (least=mΔ; oleast=ord;)
        @show ord,mΔ
    end
    oleast
end
testgh() # 1418 best order for trial function

I found out that most of the nodes on the boundary produce weights which do not really contribute anything (subnormal zeros). Thus i’m pruning them using as a threshold eps(Float64)/2.
this turned out very nice! it let’s me calculate my integrals blazing fast and with good accuracy and control

Thank you again Steven.
I’m also grateful to the author of FastGaussQuadrature.

3 Likes
4

I should also comment that passing a relative tolerance of less than machine precision to quadgk has a good chance of either never terminating in a reasonable amount of time or of hitting a DomainError if you have singular endpoints. It’s not really reasonable to demand this level of accuracy.

(It’s also suspicious to me that you are passing the same values for atol and rtol, because the two quantities have different “units”. For example, an atol=1e-16 is absurdly huge tolerance if your integrand is 1e-100 * sin(x) integrated from x=0…1, and an absurdly small tolerance if your integrand is 1e+100 * sin(x).)

2 Likes
5

Thank you again!
i agree completely, at the beggining i was testing more complicated code, that quadGK call was indeed not optimal for that trial function.
i fixed the post, it is never good to share subpar code!
the happy coincidence is that the best order remained unchanged

Thank you again Steven