FastGaussQuadrature.jl
FastGaussQuadrature.jl copied to clipboard
Factor out exponential decay in weights
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
This is important for Fun(x -> sin(x)*exp(-x), WeightedLaguerre())
.
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 !
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).
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
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.
Both work with the LinearAlgebra.jl package
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.
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
Still, the fact that it works for tridiagonal matrices is very useful and is good to know of.