From 9d16ac2b684ac6a123291359b592dc5e0a7e6929 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 10 Dec 2019 08:21:32 -0500 Subject: [PATCH] more accurate weighted gauss via chebyshev (#42) * more accurate weighted gauss via chebyshev * added analytical result for Hermite test * whoops, only evaluate cheb polynomials in (-1,1) * update README --- Project.toml | 4 +- README.md | 8 ++-- src/adapt.jl | 28 ++++++++------ src/weightedgauss.jl | 91 +++++++++++++++++++++++++++++++++++--------- test/runtests.jl | 17 ++++++--- 5 files changed, 104 insertions(+), 44 deletions(-) diff --git a/Project.toml b/Project.toml index f5fa17c..d94e7ae 100644 --- a/Project.toml +++ b/Project.toml @@ -1,16 +1,14 @@ name = "QuadGK" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.2.0" +version = "2.3.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" [compat] DataStructures = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17" julia = "0.7, 1" -Polynomials = "0.5, 0.6" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index d1f8baf..17777de 100644 --- a/README.md +++ b/README.md @@ -33,17 +33,17 @@ integration — with a little experimentation, you may be able to decide on an a you can use `x, w = gauss(N, a, b)` to find the quadrature points `x` and weights `w`, so that `sum(f.(x) .* w)` is an `N`-point approximation to `∫f(x)dx` from `a` to `b`. -For computing many integrands of similar functions with *singularities* (and/or infinite intervals), +For computing many integrands of similar functions with *singularities*, `x, w = gauss(W, N, a, b)` function allows you to pass a *weight function* `W(x)` as the first argument, so that `sum(f.(x) .* w)` is an `N`-point approximation to `∫W(x)f(x)dx` from `a` to `b`. In this way, you can put all of the singularities etcetera into `W` and precompute an accurate quadrature rule as long as the remaining `f(x)` terms are smooth. For example, ```jl using QuadGK -x, w = gauss(x -> exp(-x) / sqrt(x), 10, 0, Inf, rtol=1e-10) +x, w = gauss(x -> exp(-x) / sqrt(x), 10, 0, 20) ``` -computes the points and weights for performing `∫exp(-x)f(x)/√x dx` integrals from `0` to `∞`, so that -there is a `1/√x` singularity in the integrand at `x=0` and an implicit singularity from the `x=∞` endpoint. For example, with `f(x) = sin(x)`, the exact answer is `0.5^0.25*sin(π/8)*√π ≈ 0.5703705559915792`. Using the points and weights above with `sum(sin.(x) .* w)`, we obtain `0.5703705399484373`, which is correct to 7 digits using only 10 `f(x)` evaluations. Obtaining similar +computes the points and weights for performing `∫exp(-x)f(x)/√x dx` integrals from `0` to `20`, so that +there is a `1/√x` singularity in the integrand at `x=0` and a rapid decay for increasing `x`. For example, with `f(x) = sin(x)`, the exact answer is `0.57037055568959477…`. Using the points and weights above with `sum(sin.(x) .* w)`, we obtain `0.5703706546318231`, which is correct to 6–7 digits using only 10 `f(x)` evaluations. Obtaining similar accuracy for the same integral from `quadgk` requires nearly 300 function evaluations. However, the `gauss` function itself computes many (`2N`) numerical integrals of your weight function (multiplied by polynomials), so this is only more efficient if your `f(x)` is very expensive or if you need diff --git a/src/adapt.jl b/src/adapt.jl index 15aba5c..c235196 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -62,33 +62,35 @@ function adapt(f, segs::Vector{T}, I, E, numevals, x,w,gw,n, atol, rtol, maxeval return (I, E) end -# transform f and the endpoints s to handle infinite intervals, if any. -function handle_infinities(f, s, n, atol, rtol, maxevals, nrm) - if eltype(s) <: Real # check for infinite or semi-infinite intervals - s1 = s[1]; s2 = s[end]; inf1 = isinf(s1); inf2 = isinf(s2) +# transform f and the endpoints s to handle infinite intervals, if any, +# and pass transformed data to workfunc(f, s, tfunc) +function handle_infinities(workfunc, f, s) + s1, s2 = s[1], s[end] + if s1 isa Real && s2 isa Real # check for infinite or semi-infinite intervals + inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return do_quadgk(t -> begin t2 = t*t; den = 1 / (1 - t2); + return workfunc(t -> begin t2 = t*t; den = 1 / (1 - t2); f(t*den) * (1+t2)*den*den; end, map(x -> isinf(x) ? copysign(one(x), x) : 2x / (1+hypot(1,2x)), s), - n, atol, rtol, maxevals, nrm) + t -> t / (1 - t^2)) end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < 0 # x = s0 - t/(1-t) - return do_quadgk(t -> begin den = 1 / (1 - t); + return workfunc(t -> begin den = 1 / (1 - t); f(s0 - t*den) * den*den; end, reverse(map(x -> 1 / (1 + 1 / (s0 - x)), s)), - n, atol, rtol, maxevals, nrm) + t -> s0 - t/(1-t)) else # x = s0 + t/(1-t) - return do_quadgk(t -> begin den = 1 / (1 - t); + return workfunc(t -> begin den = 1 / (1 - t); f(s0 + t*den) * den*den; end, map(x -> 1 / (1 + 1 / (x - s0)), s), - n, atol, rtol, maxevals, nrm) + t -> s0 + t/(1-t)) end end end end - return do_quadgk(f, s, n, atol, rtol, maxevals, nrm) + return workfunc(f, s, identity) end # Gauss-Kronrod quadrature of f from a to b to c... @@ -153,4 +155,6 @@ quadgk(f, a, b, c...; kws...) = quadgk(f, a::T,b::T,c::T...; atol=nothing, rtol=nothing, maxevals=10^7, order=7, norm=norm) where {T} = - handle_infinities(f, (a, b, c...), order, atol, rtol, maxevals, norm) + handle_infinities(f, (a, b, c...)) do f, s, _ + do_quadgk(f, s, order, atol, rtol, maxevals, norm) + end diff --git a/src/weightedgauss.jl b/src/weightedgauss.jl index 8b7c333..88144ca 100644 --- a/src/weightedgauss.jl +++ b/src/weightedgauss.jl @@ -1,7 +1,42 @@ # Constructing Gaussian quadrature weights for an # arbitrary weight function and integration bounds. -using Polynomials +# for numerical stability, we apply the usual Lanczos +# Gram–Schmidt procedure to the basis {T₀,T₁,T₂,…} of +# Chebyshev polynomials on [-1,1] rather than to the +# textbook monomial basis {1,x,x²,…}. + +# evaluate Chebyshev polynomial p(x) with coefficients a[i] +# by a Clenshaw recurrence. +function chebeval(x, a) + if length(a) ≤ 2 + length(a) == 1 && return a[1] + x * zero(a[1]) + return a[1]+x*a[2] + end + bₖ = a[end-1] + 2x*a[end] + bₖ₊₁ = oftype(bₖ, a[end]) + for j = lastindex(a)-2:-1:2 + bⱼ = a[j] + 2x*bₖ - bₖ₊₁ + bₖ, bₖ₊₁ = bⱼ, bₖ + end + return a[1] + x*bₖ - bₖ₊₁ +end + +# if a[i] are coefficients of Chebyshev series, compute the coefficients xa (in-place) +# of the series multiplied by x, using recurrence xTₙ = 0.5 (Tₙ₊₁+Tₙ₋₁) for n > 0 +function chebx!(xa, a) + resize!(xa, length(a)+1) + xa .= 0 + for n = 2:lastindex(a) + c = 0.5*a[n] + xa[n-1] += c + xa[n+1] += c + end + if !isempty(a) + xa[2] += a[1] + end + return xa +end """ gauss(W, N, a, b; rtol=sqrt(eps), quad=quadgk) @@ -22,33 +57,51 @@ via the `quad` keyword argument, which should accept arguments `quad(f,a,b,rtol= similar to `quadgk`. (This is useful if your weight function has discontinuities, in which case you might want to break up the integration interval at the discontinuities.) """ -function gauss(W, N, a::Real,b::Real; rtol::Real=sqrt(eps(typeof(float(b-a)))), quad=quadgk) +gauss(W, N, a::Real,b::Real; rtol::Real=sqrt(eps(typeof(float(b-a)))), quad=quadgk) = + handle_infinities(W, (a,b)) do W, ab, tfunc + x, w = _gauss(W, N, tfunc, ab..., rtol, quad) + return (x, w) + end + +function _gauss(W, N, tfunc, a, b, rtol, quad) # Uses the Lanczos recurrence described in Trefethen & Bau, # Numerical Linear Algebra, to find the `N`-point Gaussian quadrature - # using O(N) integrals and O(N²) operations. - α = zeros(N) - β = zeros(N+1) - T = typeof(float(b-a)) - x = Poly(T[0, 1]) - q₀ = Poly(T[0]) - wint = first(quad(W, a, b, rtol=rtol)) + # using O(N) integrals and O(N²) operations, applied to Chebyshev basis: + xscale = 0.5*(b-a) # scaling from [-1,1] to (a,b) + T = typeof(xscale) + α = zeros(T, N) + β = zeros(T, N+1) + q₀ = sizehint!(T[0], N+1) # 0 polynomial + # (note that our integrals are rescaled to [-1,1], with the Jacobian + # factor |xscale| absorbed into the definition of the inner product.) + wint = first(quad(x -> W(a + (x+1)*xscale), T(-1), T(1), rtol=rtol)) + (wint isa Real && wint > 0) || + throw(ArgumentError("weight W must be real and positive")) atol = rtol*wint - q₁ = Poly(T[T(1)/sqrt(wint)]) + q₁ = sizehint!(T[1/sqrt(wint)],N+1) # normalized constant polynomial + v = copy(q₀) for n = 1:N - v = x * q₁ - q₁v = q₁ * v - α[n] = first(quad(x -> W(x) * q₁v(x), a, b, rtol=rtol, atol=atol)) + chebx!(v, q₁) # v = x * q₁ + α[n] = first(quad(T(-1), T(1), rtol=rtol, atol=atol) do x + y = a + (x+1)*xscale + W(y) * chebeval(x, q₁) * chebeval(x, v) + end) n == N && break - v -= β[n]*q₀ + α[n]*q₁ - v² = v*v - β[n+1] = sqrt(first(quad(x -> W(x) * v²(x), a, b, rtol=rtol, atol=atol))) - q₀ = q₁ - q₁ = v / β[n+1] + for j = 1:length(q₀); v[j] -= β[n]*q₀[j]; end + for j = 1:length(q₁); v[j] -= α[n]*q₁[j]; end + β[n+1] = sqrt(first(quad(T(-1), T(1), rtol=rtol, atol=atol) do x + y = a + (x+1)*xscale + W(y) * chebeval(x, v)^2 + end)) + v .*= inv(β[n+1]) + q₀,q₁,v = q₁,v,q₀ end # TODO: handle BigFloat etcetera — requires us to update eignewt() to # support nonzero diagonal entries. E = eigen(SymTridiagonal(α, β[2:N])) - return (E.values, wint .* abs2.(E.vectors[1,:])) + w = E.vectors[1,:] + w .= (abs(xscale) * wint) .* abs2.(w) + return (tfunc.(a .+ (E.values .+ 1) .* xscale), w) end diff --git a/test/runtests.jl b/test/runtests.jl index 45ec69a..b63db91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -66,9 +66,14 @@ end @test x ≈ x′ ≈ [-0.9739065285171717, -0.8650633666889845, -0.6794095682990244, -0.4333953941292472, -0.14887433898163124, 0.14887433898163124, 0.4333953941292472, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717] @test w ≈ w′ ≈ [0.06667134430868811, 0.14945134915058064, 0.2190863625159821, 0.26926671930999635, 0.29552422471475276, 0.29552422471475276, 0.26926671930999635, 0.2190863625159821, 0.14945134915058064, 0.06667134430868811] - # Gauss–Hermite quadrature - xH,wH = @inferred gauss(x->exp(-x^2), 5, -Inf, Inf, rtol=1e-10) - @test xH ≈ [-2.020182870456085632929, -0.9585724646138185071128, 0, 0.9585724646138185071128, 2.020182870456085632929] - @test wH ≈ [0.01995324205904591320774, 0.3936193231522411598285, 0.9453087204829418812257, 0.393619323152241159829, 0.01995324205904591320774] - @test sum(xH.^4 .* wH) ≈ quadgk(x -> x^4 * exp(-x^2), -Inf, Inf)[1] -end \ No newline at end of file + # infinite intervals — similar in spirit to Gauss–Hermite quadrature, + # but due to the way we use Chebyshev polynomials there is a change of variables + # in the polynomial exactness condition + xH,wH = @inferred gauss(x->exp(-x^2), 40, -Inf, Inf, rtol=1e-10) + @test sum(xH.^4 .* wH) ≈ quadgk(x -> x^4 * exp(-x^2), -Inf, Inf)[1] ≈ 3sqrt(π)/4 + + x,w = gauss(200, 2, 3) + x′,w′ = gauss(x->1, 200, 2, 3, rtol=1e-11) + @test x ≈ x′ + @test w ≈ w′ +end