# Advice on approximating Complicated integral

**URL:** <https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235>\
**Category:** General Usage\
**Created:** [April 26, 2020, 12:38pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235 "2020-04-26T12:38:20Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [April 26, 2020, 12:38pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235/1 "2020-04-26T12:38:20Z")

</div>

hi everyone i’m trying to solve a system of fixed point equations that involve integrals of this type.  
 ![immagine](https://global.discourse-cdn.com/julialang/original/3X/1/d/1d2ffbd6eef478ccf3975ee11a744c10477d6c1e.png)

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 🤣)  
thank you!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 26, 2020, 1:57pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235/2 "2020-04-26T13:57:15Z")

</div>

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](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) 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.)

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [April 26, 2020, 11:04pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235/3 "2020-04-26T23:04:11Z")

</div>

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

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 27, 2020, 12:48pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235/4 "2020-04-27T12:48:06Z")

</div>

> [@francesco.alemanno](#):
>
> `gmoment=quadgk(x->exp(-x*x/2)*F(x),-Inf,Inf,atol=1e-16,rtol=1e-16)`

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

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [April 27, 2020, 1:03pm UTC](https://discourse.julialang.org/t/advice-on-approximating-complicated-integral/38235/5 "2020-04-27T13:03:35Z")

</div>

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