Factor out exponential decay in weights

[dlfivefifty]

Suppose we want to find
$$
\int_0^\infty f(x) exp(-x) dx
$$
But we are not given $f(x)$, instead, we are given $g(x) exp(-x)$. Now we try Gauss--Laguerre as follows:

x, w= FastGaussQuadrature.gausslaguerre(150)
g = x -> sin(x)*exp(-x)
dot(w, g.(x).*exp.(x)) ≈ 0.5 # true

Great!

Uh oh:

x, w= FastGaussQuadrature.gausslaguerre(200)
g = x -> sin(x)*exp(-x)
dot(w, g.(x).*exp.(x))  |> isnan # true

If we had a weightedgausslaguerre that returned w.*exp.(x) instead of w everything would be perfect. I suspect it's just as easy to calculate these weights as standard Gauss–Laguerre.

@popsomer

[dlfivefifty]

This is important for Fun(x -> sin(x)*exp(-x), WeightedLaguerre()).

[popsomer]

This should be extremely easy to adjust in pull request #53 for such large n, as laguerreExp(n,alpha) explicitly multiplies with exp(-x) or more specifically with the weight x^alpha*exp(-x). It would only depend on which syntax would be preferred for this in gausslaguerre.jl. Special thanks to @daanhb for speeding up that code !

[dlfivefifty]

Awesome! I tried to get it working via Golub--Welsch, but it's a lot more complicated than I thought (and not obvious that it's at all possible).

[dlfivefifty]

I came up with a fairly easy workaround use BigFloat with Golub–Welsch, via https://github.com/andreasnoack/LinearAlgebra.jl:

using LinearAlgebra

function laguerreGW( n::Integer, alpha )
# Calculate Gauss-Laguerre nodes and weights based on Golub-Welsch

    alph = 2*(1:n) .+ (alpha-1)           # 3-term recurrence coeffs
    beta = sqrt.( (1:n-1).*(alpha .+ (1:n-1) ) )
    T = SymTridiagonal(Vector(alph), beta)  # Jacobi matrix
    x, V = eig( T)                  # eigenvalue decomposition
    w = gamma(alpha+1)*V[1,:].^2     # Quadrature weights
    x, vec(w)

end


let (x,w) = (Dict{Tuple{Int,BigFloat},Vector{BigFloat}}(), Dict{Tuple{Int,BigFloat},Vector{BigFloat}}())
    global function gl(n, α_in=0)
        α = BigFloat(α_in)
        haskey(x,(n,α)) && return (x[(n,α)], w[(n,α)])
        xx,ww = laguerreGW(n, α)
        x[(n,α)] = xx
        w[(n,α)] = ww
        xx,ww
    end
end

[daanhb]

Neat. I wasn't aware that eig on BigFloat tridiagonal matrices would work, since it doesn't for full matrices.

Modifying the expansions would indeed be simple and, especially for Laguerre, it makes sense to have a way to return the quadrature weights without the weight function in them. We did not deduce many terms in this regime though, since the weights rapidly underflow. The absolute error of the weights in the Airy region is extremely tiny, but the relative error may still be quite significant. We are experimenting with different implementations of the expansions over in PolynomialAsympttotics.jl.

[dlfivefifty]

Both work with the LinearAlgebra.jl package

[daanhb]

Perhaps I'm not using the right version then. I receive a MethodError related to eigfact for a dense BigFloat. It does work for tridiagonal matrices. I'll look into it.

[dlfivefifty]

Oh I didn’t actually check

On second thought it makes sense that it’s not implemented for non-Symmetric matrices since that’s a harder problem

[daanhb]

Still, the fact that it works for tridiagonal matrices is very useful and is good to know of.