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

@@ -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"

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

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

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

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
+    # 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