From 28409925aa63862447364cb93114ff7e93f303a0 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Tue, 28 Nov 2023 11:55:02 -0500 Subject: [PATCH 1/8] implement boundary asymptotics with only low order terms (tests will fail by a digit without higher order terms) --- src/gaussjacobi.jl | 225 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 208 insertions(+), 17 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index 435c8cb..358163a 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -174,7 +174,7 @@ function weightsConstant(n, α, β) end function jacobi_asy(n::Integer, α::Float64, β::Float64) - # ASY Compute nodes and weights using asymptotic formulae. + # ASY Compute nodes and weights using asymptotic formulae. # Determine switch between interior and boundary regions: nbdy = 10 @@ -185,12 +185,12 @@ function jacobi_asy(n::Integer, α::Float64, β::Float64) x, w = asy1(n, α, β, nbdy) # Boundary formula (right): - xbdy = boundary(n, α, β, nbdy) + xbdy = asy2(n, α, β, nbdy) x[bdyidx1], w[bdyidx1] = xbdy # Boundary formula (left): if α ≠ β - xbdy = boundary(n, β, α, nbdy) + xbdy = asy2(n, β, α, nbdy) end x[bdyidx2] = -xbdy[1] w[bdyidx2] = xbdy[2] @@ -350,10 +350,10 @@ function feval_asy1(n::Integer, α::Float64, β::Float64, t::AbstractVector, idx g = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880, 5246819/75246796800, -534703531/902961561600, -4483131259/86684309913600, 432261921612371/514904800886784000] - f(g,z) = dot(g, [1;cumprod(ones(9)./z)]) + f(z) = dot(g, [1;cumprod(ones(9)./z)]) # Float constant C, C2 - C = 2*p2*(f(g,n+α)*f(g,n+β)/f(g,2n+α+β))/π + C = 2*p2*(f(n+α)*f(n+β)/f(2n+α+β))/π C2 = C*(α+β+2n)*(α+β+1+2n)/(4*(α+n)*(β+n)) vals = C*S @@ -367,39 +367,230 @@ function feval_asy1(n::Integer, α::Float64, β::Float64, t::AbstractVector, idx return vals, ders end -function boundary(n::Integer, α::Float64, β::Float64, npts::Integer) -# Algorithm for computing nodes and weights near the boundary. +function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) + # Algorithm for computing nodes and weights near the boundary. # Use Newton iterations to find the first few Bessel roots: smallK = min(30, npts) jk = approx_besselroots(α, smallK) + # Use asy formula for larger ones (See NIST 10.21.19, Olver 1974 p247) + if (npts > smallK) + μ = 4*α^2 + a8 = 8*(transpose(length(jk)+1:npts)+.5*α-.25)*pi + jk2 = .125*a8-(μ-1)./a8 - 4*(μ-1)*(7*μ-31)/3 ./ a8.^3 - + 32*(μ-1)*(83*μ.^2-982*μ+3779)/15 ./ a8.^5 - + 64*(μ-1)*(6949*μ^3-153855*μ^2+1585743*μ-6277237)/105 ./ a8.^7 + jk = [jk; jk2] + end + jk = real(jk[1:npts]) # Approximate roots via asymptotic formula: (see Olver 1974) phik = jk/(n + .5*(α + β + 1)) - x = cos.( phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(8*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + 0.5*(α + β + 1))^2 ) + t = phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(8*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + 0.5*(α + β + 1))^2 # Newton iteration: for _ in 1:10 - vals, ders = innerjacobi_rec(n, x, α, β) # Evaluate via asymptotic formula. - dx = -vals./ders # Newton update. - x += dx # Next iterate. - if norm(dx,Inf) < sqrt(eps(Float64))/200 + vals, ders = feval_asy2(n, α, β, t) # Evaluate via asymptotic formula. + dt = vals./ders # Newton update. + t += dt # Next iterate. + if norm(dt,Inf) < sqrt(eps(Float64))/200 break end end - vals, ders = innerjacobi_rec(n, x, α, β) # Evaluate via asymptotic formula. - dx = -vals./ders - x += dx + vals, ders = feval_asy2(n, α, β, t) # Evaluate via asymptotic formula. + dt = vals./ders + t += dt # flip: - x = reverse(x) + t = reverse(t) ders = reverse(ders) # Revert to x-space: - w = 1 ./ ((1 .- x.^2) .* ders.^2) + x = cos.(t) + w = transpose(1 ./ ders.^2) + v = sin.(t)./ders + return x, w end +""" +Evaluate the boundary asymptotic formula at x = cos(t). +Assumption: +* `length(t) == n ÷ 2` +""" +function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector) + # Number of elements in t: + N = length(t) + + rho = n + .5*(α + β + 1) + rho2 = n + .5*(α + β - 1) + A = (.25 - α^2) + B = (.25 - β^2) + + # Evaluate the Bessel functions: + Ja = besselj.(α, rho*t) + Jb = besselj.(α + 1, rho*t) + Jbb = besselj.(α + 1, rho2*t) + Jab = besselj.(α, rho2*t) + + # Evaluate functions for recursive definition of coefficients: + gt = A*(cot.(t/2) .- (2 ./ t)) .- B*tan.(t/2) + gtdx = A*(2 ./ t.^2 .- .5*csc.(t/2).^2) .- .5*B*sec.(t/2).^2 + tB0 = .25*gt + A10 = α*(A+3*B)/24 + A1 = gtdx/8 .- (1+2*α)/8*gt./t .- gt.^2/32 .- A10 + # # Higher terms: + # tB1t = tB1(t) + # A2t = A2(t) + + vals = Ja + Jb.*tB0/rho + Ja.*A1/rho^2 # + Jb.*tB1t/rho^3 + Ja.*A2t/rho^4 + vals2 = Jab + Jbb.*tB0/rho2 + Jab.*A1/rho2^2 # + Jbb.*tB1t/rho2^3 + Jab.*A2t/rho2^4 + + # # Higher terms (not needed for n > 1000). + # tB2t = tB2(t) + # A3t = A3(t) + # vals += Jb.*tB2t/rho^5 + Ja.*A3t/rho^6 + # vals2 += Jbb.*tB2t/rho2^5 + Jab.*A3t/rho2^6 + + # Constant out the front (Computed accurately!) + ds = .5*(α^2)/n + s = ds + jj = 1 + while abs(ds/s) > eps(Float64)/10 + jj = jj+1 + ds = -(jj-1)/(jj+1)/n*(ds*α) + s = s + ds + end + p2 = exp(s)*sqrt((n+α)/n)*(n/rho)^α + g = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880, + 5246819/75246796800, -534703531/902961561600, + -4483131259/86684309913600, 432261921612371/514904800886784000] + f(z) = dot(g, [1;cumprod(ones(9)./z)]) + C = p2*(f(n+α)/f(n))/sqrt(2) + + # Scaling: + valstmp = C*vals + denom = sin.(t/2).^(α+.5).*cos.(t/2).^(β+.5) + vals = sqrt.(t).*valstmp./denom + + # Relation for derivative: + C2 = C*n/(n+α)*(rho/rho2)^α + ders = (n*(α-β .- (2*n+α+β)*cos.(t)).*valstmp .+ 2*(n+α)*(n+β)*C2*vals2)/(2*n+α+β) + ders = ders.*(sqrt.(t)./(denom.*sin.(t))) + + return vals, ders +end + +function asy2_higherterms(α, β, theta, n) + # Higher-order terms for boundary asymptotic series. + # Compute the higher order terms in asy2 boundary formula. See [2]. + + # These constants are more useful than α and β: + A = (0.25 - α^2) + B = (0.25 - β^2) + + # For now, just work on half of the domain: + c = max(max(theta), .5) + if n < 30 + N = ceil(40 - n) + elseif n >= 30 && c > pi/2-.5 + N = 15 + else + N = 10 + end + Nm1 = N - 1 + + # Scaled 2nd-kind Chebyshev points and barycentric weights: + t = .5*c*transpose(sin.(pi*(-Nm1:2:Nm1)/(2*Nm1)) + 1) + v = [.5; ones(Nm1,1)] + v[2:2:end] = -1 + v[end] = .5*v[end] + + # The g's: + g = A*(cot.(t/2) - 2. /t) - B*tan.(t/2) + gp = A*(2 ./ t.^2 - .5*csc.(t/2).^2) - .5*(.25-β^2)*sec.(t/2).^2 + gpp = A*(-4 ./ t.^3 + .25*sin.(t).*csc.(t/2).^4) - 4*B*sin.(t/2).^4 .* csc.(t).^3 + g[1] = 0 + gp[1] = -A/6-.5*B + gpp[1] = 0 + + # B0: + B0 = .25*g./t + B0p = .25*(gp./t-g./t.^2) + B0[1] = .25*(-A/6-.5*B) + B0p[1] = 0 + + # A1: + A10 = α*(A+3*B)/24 + A1 = .125*gp - (1+2*α)/2*B0 - g.^2/32 - A10 + A1p = .125*gpp - (1+2*α)/2*B0p - gp.*g/16 + A1p_t = A1p./t + A1p_t[1] = -A/720 - A^2/576 - A*B/96 - B^2/64 - B/48 + α*(A/720 + B/48) + + # Make f accurately: (Taylor series approx for small t) + fcos = B./(2*cos.(t/2)).^2 + f = -A*(1/12 + t.^2/240+t.^4/6048 + t.^6/172800 + t.^8/5322240 + + 691*t.^10/118879488000 + t.^12/5748019200 + + 3617*t.^14/711374856192000 + 43867*t.^16/300534953951232000) + idx = t > .5 + ti = t[idx] + f[idx] = A.*(1 ./ ti.^2 - 1 ./ (2*sin(ti/2)).^2) + f = f - fcos + + # Integrals for B1: (Note that N isn't large, so we don't need to be fancy). + C = chebcolloc2.cumsummat(N)*(.5*c) + D = chebcolloc2.diffmat(N)*(2/c) + I = (C*A1p_t) + J = (C*(f.*A1)) + + # B1: + tB1 = -.5*A1p - (.5+α)*I + .5*J + tB1[1] = 0 + B1 = tB1./t + B1[1] = A/720 + A^2/576 + A*B/96 + B^2/64 + B/48 + + α*(A^2/576 + B^2/64 + A*B/96) - α^2*(A/720 + B/48) + + # A2: + K = C*(f.*tB1) + A2 = .5*(D*tB1) - (.5+α)*B1 - .5*K + A2 = A2 - A2[1] + + # if ( nargout < 3 ) + # # Make function for output + # tB1(theta) = bary(theta,tB1,t,v) + # A2(theta) = bary(theta,A2,t,v) + # return + # end + + # A2p: + A2p = D*A2 + A2p = A2p - A2p[1] + A2p_t = A2p./t + # Extrapolate point at t = 0: + w = pi/2-t[2:end] + w[2:2:end] = -w[2:2:end] + w[end] = .5*w[end] + A2p_t[1] = sum(w.*A2p_t[2:end])/sum(w) + + # B2: + tB2 = -.5*A2p - (.5+α)*(C*A2p_t) + .5*C*(f.*A2) + B2 = tB2./t + # Extrapolate point at t = 0: + B2[1] = sum(w.*B2[2:end])/sum(w) + + # A3: + K = C*(f.*tB2) + A3 = .5*(D*tB2) - (.5+α)*B2 - .5*K + A3 = A3 - A3[1] + + tB1(theta) = bary(theta, tB1, t, v) + A2(theta) = bary(theta, A2, t, v) + tB2(theta) = bary(theta, tB2, t, v) + A3(theta) = bary(theta, A3, t, v) + + return [tB1, A2, tB2, A3, tB3, A4] +end + function jacobi_jacobimatrix(n, α, β) s = α + β ii = 2:n-1 From d140b63b69ee7c61a1a2c513e4d1fe911969e036 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Thu, 30 Nov 2023 15:58:13 -0500 Subject: [PATCH 2/8] use saved chebfun integration and differentiation matrices --- src/gaussjacobi.jl | 188 ++++++++++++++++++++++++++------------------- 1 file changed, 111 insertions(+), 77 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index 358163a..6458b70 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -264,9 +264,6 @@ function feval_asy1(n::Integer, α::Float64, β::Float64, t::AbstractVector, idx # Number of elements in t: N = length(t) - # Some often used vectors/matrices: - onesM = ones(M) - # The sine and cosine terms: A = repeat((2n+α+β).+(1:M),1,N).*repeat(t',M)/2 .- (α+1/2)*π/2 # M × N matrix cosA = cos.(A) @@ -386,18 +383,21 @@ function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) # Approximate roots via asymptotic formula: (see Olver 1974) phik = jk/(n + .5*(α + β + 1)) - t = phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(8*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + 0.5*(α + β + 1))^2 + t = phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(2*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + .5*(α + β + 1))^2 + + # Compute higher terms: + higherterms = asy2_higherterms(α, β, t, n) # Newton iteration: for _ in 1:10 - vals, ders = feval_asy2(n, α, β, t) # Evaluate via asymptotic formula. + vals, ders = feval_asy2(n, α, β, t, higherterms) # Evaluate via asymptotic formula. dt = vals./ders # Newton update. t += dt # Next iterate. if norm(dt,Inf) < sqrt(eps(Float64))/200 break end end - vals, ders = feval_asy2(n, α, β, t) # Evaluate via asymptotic formula. + vals, ders = feval_asy2(n, α, β, t, higherterms) # Evaluate via asymptotic formula. dt = vals./ders t += dt @@ -418,11 +418,8 @@ Evaluate the boundary asymptotic formula at x = cos(t). Assumption: * `length(t) == n ÷ 2` """ -function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector) - # Number of elements in t: - N = length(t) - - rho = n + .5*(α + β + 1) +function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, higherterms) + rho = n + .5*(α + β + 1) rho2 = n + .5*(α + β - 1) A = (.25 - α^2) B = (.25 - β^2) @@ -439,18 +436,21 @@ function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector) tB0 = .25*gt A10 = α*(A+3*B)/24 A1 = gtdx/8 .- (1+2*α)/8*gt./t .- gt.^2/32 .- A10 - # # Higher terms: - # tB1t = tB1(t) - # A2t = A2(t) - - vals = Ja + Jb.*tB0/rho + Ja.*A1/rho^2 # + Jb.*tB1t/rho^3 + Ja.*A2t/rho^4 - vals2 = Jab + Jbb.*tB0/rho2 + Jab.*A1/rho2^2 # + Jbb.*tB1t/rho2^3 + Jab.*A2t/rho2^4 + # Higher terms: + tB1, A2, tB2, A3 = higherterms + tB1t = tB1(t) + A2t = A2(t) + + # VALS: + vals = Ja + Jb.*tB0/rho + Ja.*A1/rho^2 + Jb.*tB1t/rho^3 + Ja.*A2t/rho^4 + # DERS: + vals2 = Jab + Jbb.*tB0/rho2 + Jab.*A1/rho2^2 + Jbb.*tB1t/rho2^3 + Jab.*A2t/rho2^4 - # # Higher terms (not needed for n > 1000). - # tB2t = tB2(t) - # A3t = A3(t) - # vals += Jb.*tB2t/rho^5 + Ja.*A3t/rho^6 - # vals2 += Jbb.*tB2t/rho2^5 + Jab.*A3t/rho2^6 + # Higher terms (not needed for n > 1000). + tB2t = tB2(t) + A3t = A3(t) + vals .+= Jb.*tB2t/rho^5 + Ja.*A3t/rho^6 + vals2 .+= Jbb.*tB2t/rho2^5 + Jab.*A3t/rho2^6 # Constant out the front (Computed accurately!) ds = .5*(α^2)/n @@ -475,13 +475,13 @@ function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector) # Relation for derivative: C2 = C*n/(n+α)*(rho/rho2)^α - ders = (n*(α-β .- (2*n+α+β)*cos.(t)).*valstmp .+ 2*(n+α)*(n+β)*C2*vals2)/(2*n+α+β) + ders = (n*(α-β .- (2n+α+β)*cos.(t)).*valstmp .+ 2*(n+α)*(n+β)*C2*vals2)/(2n+α+β) ders = ders.*(sqrt.(t)./(denom.*sin.(t))) return vals, ders end -function asy2_higherterms(α, β, theta, n) +function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector, n::Integer) # Higher-order terms for boundary asymptotic series. # Compute the higher order terms in asy2 boundary formula. See [2]. @@ -490,105 +490,89 @@ function asy2_higherterms(α, β, theta, n) B = (0.25 - β^2) # For now, just work on half of the domain: - c = max(max(theta), .5) - if n < 30 - N = ceil(40 - n) - elseif n >= 30 && c > pi/2-.5 - N = 15 - else - N = 10 - end - Nm1 = N - 1 + c = max(maximum(theta), .5) # Scaled 2nd-kind Chebyshev points and barycentric weights: - t = .5*c*transpose(sin.(pi*(-Nm1:2:Nm1)/(2*Nm1)) + 1) - v = [.5; ones(Nm1,1)] - v[2:2:end] = -1 - v[end] = .5*v[end] + t = .5*c*(sin.(pi*(-(N-1):2:(N-1))/(2*(N-1))) .+ 1) + v = [.5; ones(N-1,1)] + v[2:2:end] .= -1 + v[end] *= .5 # The g's: - g = A*(cot.(t/2) - 2. /t) - B*tan.(t/2) - gp = A*(2 ./ t.^2 - .5*csc.(t/2).^2) - .5*(.25-β^2)*sec.(t/2).^2 + g = A*(cot.(t/2) - 2 ./t) - B*tan.(t/2) + gp = A*(2 ./ t.^2 - .5*csc.(t/2).^2) - .5*(.25-β^2)*sec.(t/2).^2 gpp = A*(-4 ./ t.^3 + .25*sin.(t).*csc.(t/2).^4) - 4*B*sin.(t/2).^4 .* csc.(t).^3 - g[1] = 0 - gp[1] = -A/6-.5*B + g[1] = 0 + gp[1] = -A/6-.5*B gpp[1] = 0 # B0: - B0 = .25*g./t - B0p = .25*(gp./t-g./t.^2) - B0[1] = .25*(-A/6-.5*B) + B0 = .25*g./t + B0p = .25*(gp./t - g./t.^2) + B0[1] = .25*(-A/6-.5*B) B0p[1] = 0 # A1: - A10 = α*(A+3*B)/24 - A1 = .125*gp - (1+2*α)/2*B0 - g.^2/32 - A10 - A1p = .125*gpp - (1+2*α)/2*B0p - gp.*g/16 + A10 = α*(A+3*B)/24 + A1 = .125*gp .- (1+2*α)/2*B0 .- g.^2/32 .- A10 + A1p = .125*gpp .- (1+2*α)/2*B0p .- gp.*g/16 A1p_t = A1p./t A1p_t[1] = -A/720 - A^2/576 - A*B/96 - B^2/64 - B/48 + α*(A/720 + B/48) # Make f accurately: (Taylor series approx for small t) fcos = B./(2*cos.(t/2)).^2 - f = -A*(1/12 + t.^2/240+t.^4/6048 + t.^6/172800 + t.^8/5322240 + + f = -A*(1/12 .+ t.^2/240+t.^4/6048 + t.^6/172800 + t.^8/5322240 + 691*t.^10/118879488000 + t.^12/5748019200 + 3617*t.^14/711374856192000 + 43867*t.^16/300534953951232000) - idx = t > .5 - ti = t[idx] - f[idx] = A.*(1 ./ ti.^2 - 1 ./ (2*sin(ti/2)).^2) - f = f - fcos + idx = t .> .5 + f[idx] = A.*(1 ./ t[idx].^2 - 1 ./ (2*sin.(t[idx]/2)).^2) + f .-= fcos # Integrals for B1: (Note that N isn't large, so we don't need to be fancy). - C = chebcolloc2.cumsummat(N)*(.5*c) - D = chebcolloc2.diffmat(N)*(2/c) + C = (.5*c)*cumsummatN + D = (2/c)*diffmatN I = (C*A1p_t) J = (C*(f.*A1)) # B1: - tB1 = -.5*A1p - (.5+α)*I + .5*J + tB1 = -.5*A1p - (.5+α)*I + .5*J tB1[1] = 0 - B1 = tB1./t - B1[1] = A/720 + A^2/576 + A*B/96 + B^2/64 + B/48 + + B1 = tB1./t + B1[1] = A/720 + A^2/576 + A*B/96 + B^2/64 + B/48 + α*(A^2/576 + B^2/64 + A*B/96) - α^2*(A/720 + B/48) # A2: - K = C*(f.*tB1) + K = C*(f.*tB1) A2 = .5*(D*tB1) - (.5+α)*B1 - .5*K - A2 = A2 - A2[1] - - # if ( nargout < 3 ) - # # Make function for output - # tB1(theta) = bary(theta,tB1,t,v) - # A2(theta) = bary(theta,A2,t,v) - # return - # end + A2 .-= A2[1] # A2p: A2p = D*A2 - A2p = A2p - A2p[1] + A2p .-= A2p[1] A2p_t = A2p./t # Extrapolate point at t = 0: - w = pi/2-t[2:end] + w = pi/2 .- t[2:end] w[2:2:end] = -w[2:2:end] - w[end] = .5*w[end] + w[end] = .5*w[end] A2p_t[1] = sum(w.*A2p_t[2:end])/sum(w) # B2: tB2 = -.5*A2p - (.5+α)*(C*A2p_t) + .5*C*(f.*A2) - B2 = tB2./t + B2 = tB2./t # Extrapolate point at t = 0: B2[1] = sum(w.*B2[2:end])/sum(w) # A3: - K = C*(f.*tB2) + K = C*(f.*tB2) A3 = .5*(D*tB2) - (.5+α)*B2 - .5*K - A3 = A3 - A3[1] + A3 .-= A3[1] - tB1(theta) = bary(theta, tB1, t, v) - A2(theta) = bary(theta, A2, t, v) - tB2(theta) = bary(theta, tB2, t, v) - A3(theta) = bary(theta, A3, t, v) + tB1f(theta) = bary(theta, tB1, t, v) + A2f(theta) = bary(theta, A2, t, v) + tB2f(theta) = bary(theta, tB2, t, v) + A3f(theta) = bary(theta, A3, t, v) - return [tB1, A2, tB2, A3, tB3, A4] + return (tB1f, A2f, tB2f, A3f) end function jacobi_jacobimatrix(n, α, β) @@ -617,3 +601,53 @@ function jacobi_gw(n::Integer, α, β) w = V[1,:].^2 .* jacobimoment(α, β) return x, w end + +function bary(x, fvals, xk, vk) + # barycentric interpolation taken from chebfun/bary.m + + # Initialise return value: + fx = zeros(length(x)) + + # Loop: + for j in eachindex(x) + xx = vk ./ (x[j] .- xk) + fx[j] = dot(xx, fvals) / sum(xx) + end + + # Try to clean up NaNs: + for k in findall(isnan.(fx)) + index = findfirst(x[k] .== xk) # Find the corresponding node + if !isempty(index) + fx[k] = fvals[index] # Set entries to the correct value + end + end + + return fx +end + +# Chebyshev type 2 integration and differentiation matrices from chebfun +const N = 10 +const cumsummatN = [ + 0 0 0 0 0 0 0 0 0 0; + 0.019080722834519 0.0496969890549313 -0.0150585059796021 0.0126377679164575 -0.0118760811432484 0.0115424841953298 -0.0113725236133433 0.0112812076497144 -0.011235316890839 0.00561063519017238; + 0.000812345683614654 0.14586999854807 0.0976007154946748 -0.0146972757610091 0.00680984376276729 -0.00401953146146086 0.00271970678005437 -0.00205195604894289 0.00172405556686793 -0.000812345683614662; + 0.017554012345679 0.103818185816131 0.249384588781868 0.149559082892416 -0.0321899366961563 0.0210262631520163 -0.0171075837742504 0.0153341224604243 -0.0145160806571407 0.00713734567901234; + 0.00286927716087872 0.136593368810421 0.201074970443365 0.339479954969535 0.164397864607267 -0.0260484364615523 0.0127399306249393 -0.00815620454308202 0.00627037388217603 -0.00286927716087872; + 0.0152149561732244 0.110297082689861 0.233440527881186 0.289200104648429 0.369910942265696 0.179464641196877 -0.0375399196961666 0.0242093528947391 -0.0200259122383839 0.00947640185146695; + 0.00520833333333334 0.131083537229178 0.20995020087768 0.319047619047619 0.322836242652128 0.376052442500301 0.152380952380952 -0.024100265443764 0.0127492707559062 -0.00520833333333333; + 0.0131580246959603 0.114843401005169 0.227336279387047 0.299220328493314 0.347882037265605 0.337052662041377 0.316637311034378 0.12768360784343 -0.0293025419760333 0.011533333328731; + 0.00673504382217329 0.127802773462876 0.21400311568839 0.313312558886712 0.332320021608814 0.355738586947393 0.289302267356911 0.240342829317707 0.0668704675171058 -0.00673504382217329; + 0.0123456790123457 0.116567456572037 0.225284323338104 0.301940035273369 0.343862505804144 0.343862505804144 0.301940035273369 0.225284323338104 0.116567456572037 0.0123456790123457 + ] +const diffmatN = [ + -27.1666666666667 33.1634374775264 -8.54863217041303 4 -2.42027662546121 1.70408819104185 -1.33333333333333 1.13247433143179 -1.03109120412576 0.5; + -8.29085936938159 4.01654328417507 5.75877048314363 -2.27431608520652 1.30540728933228 -0.898197570222574 0.694592710667722 -0.586256827714545 0.532088886237956 -0.257772801031441; + 2.13715804260326 -5.75877048314363 0.927019729872654 3.75877048314364 -1.68805925749197 1.06417777247591 -0.789861687269397 0.652703644666139 -0.586256827714545 0.283118582857949; + -1 2.27431608520652 -3.75877048314364 0.333333333333335 3.06417777247591 -1.48445439793712 1 -0.789861687269397 0.694592710667722 -0.333333333333333; + 0.605069156365302 -1.30540728933228 1.68805925749197 -3.06417777247591 0.0895235543024196 2.87938524157182 -1.48445439793712 1.06417777247591 -0.898197570222574 0.426022047760462; + -0.426022047760462 0.898197570222574 -1.06417777247591 1.48445439793712 -2.87938524157182 -0.0895235543024196 3.06417777247591 -1.68805925749197 1.30540728933228 -0.605069156365302; + 0.333333333333333 -0.694592710667722 0.789861687269397 -1 1.48445439793712 -3.06417777247591 -0.333333333333335 3.75877048314364 -2.27431608520652 1; + -0.283118582857949 0.586256827714545 -0.652703644666139 0.789861687269397 -1.06417777247591 1.68805925749197 -3.75877048314364 -0.927019729872654 5.75877048314363 -2.13715804260326; + 0.257772801031441 -0.532088886237956 0.586256827714545 -0.694592710667722 0.898197570222574 -1.30540728933228 2.27431608520652 -5.75877048314363 -4.01654328417507 8.29085936938159; + -0.5 1.03109120412576 -1.13247433143179 1.33333333333333 -1.70408819104185 2.42027662546121 -4 8.54863217041303 -33.1634374775264 27.1666666666667 + ] \ No newline at end of file From f5a58373c285382cf6f57ab91b0477dd55147bad Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Mon, 4 Dec 2023 16:04:26 -0500 Subject: [PATCH 3/8] 15 digit accuracy thanks to besselTaylor --- src/gaussjacobi.jl | 99 +++++++++++++++++++++++++++++----------- test/test_gaussjacobi.jl | 9 ++++ 2 files changed, 81 insertions(+), 27 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index 6458b70..ab136f1 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -174,7 +174,7 @@ function weightsConstant(n, α, β) end function jacobi_asy(n::Integer, α::Float64, β::Float64) - # ASY Compute nodes and weights using asymptotic formulae. + # ASY Compute nodes and weights using asymptotic formulae. # Determine switch between interior and boundary regions: nbdy = 10 @@ -381,12 +381,12 @@ function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) end jk = real(jk[1:npts]) - # Approximate roots via asymptotic formula: (see Olver 1974) + # Approximate roots via asymptotic formula: (see Olver 1974, NIST, 18.16.8) phik = jk/(n + .5*(α + β + 1)) - t = phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(2*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + .5*(α + β + 1))^2 + t = phik .+ ((α^2-0.25).*(1 .- phik.*cot.(phik))./(2*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + .5*(α + β + 1))^2 # Compute higher terms: - higherterms = asy2_higherterms(α, β, t, n) + higherterms = asy2_higherterms(α, β, t) # Newton iteration: for _ in 1:10 @@ -402,13 +402,12 @@ function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) t += dt # flip: - t = reverse(t) + t = reverse(t) ders = reverse(ders) # Revert to x-space: x = cos.(t) - w = transpose(1 ./ ders.^2) - v = sin.(t)./ders + w = 1 ./ ders.^2 return x, w end @@ -418,7 +417,7 @@ Evaluate the boundary asymptotic formula at x = cos(t). Assumption: * `length(t) == n ÷ 2` """ -function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, higherterms) +function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, higherterms::HT) where HT <: Tuple{<:Function, <:Function, <:Function, <:Function} rho = n + .5*(α + β + 1) rho2 = n + .5*(α + β - 1) A = (.25 - α^2) @@ -428,14 +427,14 @@ function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, hig Ja = besselj.(α, rho*t) Jb = besselj.(α + 1, rho*t) Jbb = besselj.(α + 1, rho2*t) - Jab = besselj.(α, rho2*t) + Jab = besselTaylor(-t, rho*t, α) # Evaluate functions for recursive definition of coefficients: gt = A*(cot.(t/2) .- (2 ./ t)) .- B*tan.(t/2) gtdx = A*(2 ./ t.^2 .- .5*csc.(t/2).^2) .- .5*B*sec.(t/2).^2 tB0 = .25*gt A10 = α*(A+3*B)/24 - A1 = gtdx/8 .- (1+2*α)/8*gt./t .- gt.^2/32 .- A10 + A1 = gtdx/8 .- (1+2*α)/8*gt./t .- gt.^2/32 .- A10 # Higher terms: tB1, A2, tB2, A3 = higherterms tB1t = tB1(t) @@ -470,18 +469,18 @@ function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, hig # Scaling: valstmp = C*vals - denom = sin.(t/2).^(α+.5).*cos.(t/2).^(β+.5) - vals = sqrt.(t).*valstmp./denom + denom = sin.(t/2).^(α+.5).*cos.(t/2).^(β+.5) + vals = sqrt.(t).*valstmp./denom # Relation for derivative: - C2 = C*n/(n+α)*(rho/rho2)^α - ders = (n*(α-β .- (2n+α+β)*cos.(t)).*valstmp .+ 2*(n+α)*(n+β)*C2*vals2)/(2n+α+β) - ders = ders.*(sqrt.(t)./(denom.*sin.(t))) + C2 = C*n/(n+α)*(rho/rho2)^α + ders = (n*(α-β .- (2n+α+β)*cos.(t)).*valstmp .+ 2*(n+α)*(n+β)*C2*vals2)/(2n+α+β) + ders .*= sqrt.(t)./(denom.*sin.(t)) return vals, ders end -function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector, n::Integer) +function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector) # Higher-order terms for boundary asymptotic series. # Compute the higher order terms in asy2 boundary formula. See [2]. @@ -493,8 +492,8 @@ function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector, n::In c = max(maximum(theta), .5) # Scaled 2nd-kind Chebyshev points and barycentric weights: - t = .5*c*(sin.(pi*(-(N-1):2:(N-1))/(2*(N-1))) .+ 1) - v = [.5; ones(N-1,1)] + t = .5*c*(sin.(pi*(-(nc-1):2:(nc-1))/(2*(nc-1))) .+ 1) + v = [.5; ones(nc-1)] v[2:2:end] .= -1 v[end] *= .5 @@ -529,8 +528,8 @@ function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector, n::In f .-= fcos # Integrals for B1: (Note that N isn't large, so we don't need to be fancy). - C = (.5*c)*cumsummatN - D = (2/c)*diffmatN + C = (.5*c)*cumsummat + D = (2/c)*diffmat I = (C*A1p_t) J = (C*(f.*A1)) @@ -572,7 +571,7 @@ function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector, n::In tB2f(theta) = bary(theta, tB2, t, v) A3f(theta) = bary(theta, A3, t, v) - return (tB1f, A2f, tB2f, A3f) + return tB1f, A2f, tB2f, A3f end function jacobi_jacobimatrix(n, α, β) @@ -602,8 +601,8 @@ function jacobi_gw(n::Integer, α, β) return x, w end -function bary(x, fvals, xk, vk) - # barycentric interpolation taken from chebfun/bary.m +function bary(x::AbstractVector, fvals::AbstractVector, xk::AbstractVector, vk::AbstractVector) + # simple barycentric interpolation routine adapted from chebfun/bary.m # Initialise return value: fx = zeros(length(x)) @@ -625,9 +624,55 @@ function bary(x, fvals, xk, vk) return fx end -# Chebyshev type 2 integration and differentiation matrices from chebfun -const N = 10 -const cumsummatN = [ +function besselTaylor(t::AbstractVector, z::AbstractVector, α::Float64) + # Accurate evaluation of Bessel function J_A for asy2. (See [2].) + # evaluates J_A(Z+T) by a Taylor series expansion about Z. + + npts = length(t) + kmax = min(ceil(Int64, abs(log(eps(Float64))/log(norm(t, Inf)))), 30) + H = t' .^ (0:kmax) + # Compute coeffs in Taylor expansions about z (See NIST 10.6.7) + nu = ones(Int64, length(z)) * (-kmax:kmax)' + JK = z * ones(2kmax+1)' + Bjk = besselj.(α .+ nu, JK) + nck = nck_mat(floor(Int64, 1.25*kmax)) # nchoosek + AA = [Bjk[:,kmax+1] zeros(npts, kmax)] + fact = 1 + for k = 1:kmax + sgn = 1 + for l = 0:k + AA[:,k+1] = AA[:,k+1] + sgn*nck[k,l+1]*Bjk[:,kmax+2*l-k+1] + sgn = -sgn + end + fact = k*fact + AA[:,k+1] ./= 2^k * fact + end + # Evaluate Taylor series: + Ja = zeros(npts) + for k = 1:npts + Ja[k] = dot(AA[k,:], H[:,k]) + end + + return Ja +end + +function nck_mat(n::Integer) + # almost triangular matrix storing n choose k + M = zeros(Int64, n-1, n) + M[:,1] .= 1 + M[1,2] = 1 + for i=2:n-1 + for j=2:i+1 + M[i,j] = M[i-1,j-1] + M[i-1,j] + end + end + return M +end + +# 10-point Chebyshev type 2 integration and differentiation matrices +# computed using Chebfun +const nc = 10 +const cumsummat = [ 0 0 0 0 0 0 0 0 0 0; 0.019080722834519 0.0496969890549313 -0.0150585059796021 0.0126377679164575 -0.0118760811432484 0.0115424841953298 -0.0113725236133433 0.0112812076497144 -0.011235316890839 0.00561063519017238; 0.000812345683614654 0.14586999854807 0.0976007154946748 -0.0146972757610091 0.00680984376276729 -0.00401953146146086 0.00271970678005437 -0.00205195604894289 0.00172405556686793 -0.000812345683614662; @@ -639,7 +684,7 @@ const cumsummatN = [ 0.00673504382217329 0.127802773462876 0.21400311568839 0.313312558886712 0.332320021608814 0.355738586947393 0.289302267356911 0.240342829317707 0.0668704675171058 -0.00673504382217329; 0.0123456790123457 0.116567456572037 0.225284323338104 0.301940035273369 0.343862505804144 0.343862505804144 0.301940035273369 0.225284323338104 0.116567456572037 0.0123456790123457 ] -const diffmatN = [ +const diffmat = [ -27.1666666666667 33.1634374775264 -8.54863217041303 4 -2.42027662546121 1.70408819104185 -1.33333333333333 1.13247433143179 -1.03109120412576 0.5; -8.29085936938159 4.01654328417507 5.75877048314363 -2.27431608520652 1.30540728933228 -0.898197570222574 0.694592710667722 -0.586256827714545 0.532088886237956 -0.257772801031441; 2.13715804260326 -5.75877048314363 0.927019729872654 3.75877048314364 -1.68805925749197 1.06417777247591 -0.789861687269397 0.652703644666139 -0.586256827714545 0.283118582857949; diff --git a/test/test_gaussjacobi.jl b/test/test_gaussjacobi.jl index 2cccf94..6a5ce63 100644 --- a/test/test_gaussjacobi.jl +++ b/test/test_gaussjacobi.jl @@ -230,4 +230,13 @@ @test abs(x[7] - big(4146701176053244724446765533149623814305606864483746748725151214196203624388371)/big(10)^79) < 1e-70 @test abs(w[3] - big(2482452398859023537636458227096017466104571146686194593251519240011574719667464)/big(10)^79) < 1e-70 end + + @testset "boundary asymptotics were previously not implemented (#58,#130)" begin + _, w = gaussjacobi(2^10, 0.25, 0) + @test abs(w[end] - 3.607554904604311e-07) / 3.607554904604311e-07 < tol_w + + _, w = gaussjacobi(2^16, -0.9, 0) + @test abs(w[end] - 1.223027243345722) < tol_w + @test abs(sum(w) - 2^(-0.9+1)/(-0.9+1)) < tol + end end \ No newline at end of file From b828ffa1e78c9c6b0063333d14cf8d1c6db7fc94 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Tue, 5 Dec 2023 11:33:36 -0500 Subject: [PATCH 4/8] move constants into constants.jl --- src/constants.jl | 35 +++++++++++++++++++++++++++++++++++ src/gaussjacobi.jl | 42 +++++++++--------------------------------- 2 files changed, 44 insertions(+), 33 deletions(-) diff --git a/src/constants.jl b/src/constants.jl index a0d9177..9609e23 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -125,3 +125,38 @@ const AIRY_ROOTS = @SVector [ -12.828776752865757, -13.69148903521072, ] + +@doc raw""" +10-point Chebyshev type 2 integration matrix computed using Chebfun. + +Used for numerical integration in boundary asymptotics for Gauss-Jacobi. +""" +const CUMSUMMAT_10 = [ + 0 0 0 0 0 0 0 0 0 0; + 0.019080722834519 0.0496969890549313 -0.0150585059796021 0.0126377679164575 -0.0118760811432484 0.0115424841953298 -0.0113725236133433 0.0112812076497144 -0.011235316890839 0.00561063519017238; + 0.000812345683614654 0.14586999854807 0.0976007154946748 -0.0146972757610091 0.00680984376276729 -0.00401953146146086 0.00271970678005437 -0.00205195604894289 0.00172405556686793 -0.000812345683614662; + 0.017554012345679 0.103818185816131 0.249384588781868 0.149559082892416 -0.0321899366961563 0.0210262631520163 -0.0171075837742504 0.0153341224604243 -0.0145160806571407 0.00713734567901234; + 0.00286927716087872 0.136593368810421 0.201074970443365 0.339479954969535 0.164397864607267 -0.0260484364615523 0.0127399306249393 -0.00815620454308202 0.00627037388217603 -0.00286927716087872; + 0.0152149561732244 0.110297082689861 0.233440527881186 0.289200104648429 0.369910942265696 0.179464641196877 -0.0375399196961666 0.0242093528947391 -0.0200259122383839 0.00947640185146695; + 0.00520833333333334 0.131083537229178 0.20995020087768 0.319047619047619 0.322836242652128 0.376052442500301 0.152380952380952 -0.024100265443764 0.0127492707559062 -0.00520833333333333; + 0.0131580246959603 0.114843401005169 0.227336279387047 0.299220328493314 0.347882037265605 0.337052662041377 0.316637311034378 0.12768360784343 -0.0293025419760333 0.011533333328731; + 0.00673504382217329 0.127802773462876 0.21400311568839 0.313312558886712 0.332320021608814 0.355738586947393 0.289302267356911 0.240342829317707 0.0668704675171058 -0.00673504382217329; + 0.0123456790123457 0.116567456572037 0.225284323338104 0.301940035273369 0.343862505804144 0.343862505804144 0.301940035273369 0.225284323338104 0.116567456572037 0.0123456790123457 + ] +@doc raw""" +10-point Chebyshev type 2 differentiation matrix computed using Chebfun. + +Used for numerical differentiation in boundary asymptotics for Gauss-Jacobi. +""" +const DIFFMAT_10 = [ + -27.1666666666667 33.1634374775264 -8.54863217041303 4 -2.42027662546121 1.70408819104185 -1.33333333333333 1.13247433143179 -1.03109120412576 0.5; + -8.29085936938159 4.01654328417507 5.75877048314363 -2.27431608520652 1.30540728933228 -0.898197570222574 0.694592710667722 -0.586256827714545 0.532088886237956 -0.257772801031441; + 2.13715804260326 -5.75877048314363 0.927019729872654 3.75877048314364 -1.68805925749197 1.06417777247591 -0.789861687269397 0.652703644666139 -0.586256827714545 0.283118582857949; + -1 2.27431608520652 -3.75877048314364 0.333333333333335 3.06417777247591 -1.48445439793712 1 -0.789861687269397 0.694592710667722 -0.333333333333333; + 0.605069156365302 -1.30540728933228 1.68805925749197 -3.06417777247591 0.0895235543024196 2.87938524157182 -1.48445439793712 1.06417777247591 -0.898197570222574 0.426022047760462; + -0.426022047760462 0.898197570222574 -1.06417777247591 1.48445439793712 -2.87938524157182 -0.0895235543024196 3.06417777247591 -1.68805925749197 1.30540728933228 -0.605069156365302; + 0.333333333333333 -0.694592710667722 0.789861687269397 -1 1.48445439793712 -3.06417777247591 -0.333333333333335 3.75877048314364 -2.27431608520652 1; + -0.283118582857949 0.586256827714545 -0.652703644666139 0.789861687269397 -1.06417777247591 1.68805925749197 -3.75877048314364 -0.927019729872654 5.75877048314363 -2.13715804260326; + 0.257772801031441 -0.532088886237956 0.586256827714545 -0.694592710667722 0.898197570222574 -1.30540728933228 2.27431608520652 -5.75877048314363 -4.01654328417507 8.29085936938159; + -0.5 1.03109120412576 -1.13247433143179 1.33333333333333 -1.70408819104185 2.42027662546121 -4 8.54863217041303 -33.1634374775264 27.1666666666667 + ] \ No newline at end of file diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index ab136f1..a7b6adb 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -491,9 +491,13 @@ function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector) # For now, just work on half of the domain: c = max(maximum(theta), .5) + # Use 10 Chebyshev modes in order to use precomputed integration and + # differentiation matrices + N = 10 + # Scaled 2nd-kind Chebyshev points and barycentric weights: - t = .5*c*(sin.(pi*(-(nc-1):2:(nc-1))/(2*(nc-1))) .+ 1) - v = [.5; ones(nc-1)] + t = .5*c*(sin.(pi*(-(N-1):2:(N-1))/(2*(N-1))) .+ 1) + v = [.5; ones(N-1)] v[2:2:end] .= -1 v[end] *= .5 @@ -528,8 +532,8 @@ function asy2_higherterms(α::Float64, β::Float64, theta::AbstractVector) f .-= fcos # Integrals for B1: (Note that N isn't large, so we don't need to be fancy). - C = (.5*c)*cumsummat - D = (2/c)*diffmat + C = (.5*c)*CUMSUMMAT_10 + D = (2/c)*DIFFMAT_10 I = (C*A1p_t) J = (C*(f.*A1)) @@ -667,32 +671,4 @@ function nck_mat(n::Integer) end end return M -end - -# 10-point Chebyshev type 2 integration and differentiation matrices -# computed using Chebfun -const nc = 10 -const cumsummat = [ - 0 0 0 0 0 0 0 0 0 0; - 0.019080722834519 0.0496969890549313 -0.0150585059796021 0.0126377679164575 -0.0118760811432484 0.0115424841953298 -0.0113725236133433 0.0112812076497144 -0.011235316890839 0.00561063519017238; - 0.000812345683614654 0.14586999854807 0.0976007154946748 -0.0146972757610091 0.00680984376276729 -0.00401953146146086 0.00271970678005437 -0.00205195604894289 0.00172405556686793 -0.000812345683614662; - 0.017554012345679 0.103818185816131 0.249384588781868 0.149559082892416 -0.0321899366961563 0.0210262631520163 -0.0171075837742504 0.0153341224604243 -0.0145160806571407 0.00713734567901234; - 0.00286927716087872 0.136593368810421 0.201074970443365 0.339479954969535 0.164397864607267 -0.0260484364615523 0.0127399306249393 -0.00815620454308202 0.00627037388217603 -0.00286927716087872; - 0.0152149561732244 0.110297082689861 0.233440527881186 0.289200104648429 0.369910942265696 0.179464641196877 -0.0375399196961666 0.0242093528947391 -0.0200259122383839 0.00947640185146695; - 0.00520833333333334 0.131083537229178 0.20995020087768 0.319047619047619 0.322836242652128 0.376052442500301 0.152380952380952 -0.024100265443764 0.0127492707559062 -0.00520833333333333; - 0.0131580246959603 0.114843401005169 0.227336279387047 0.299220328493314 0.347882037265605 0.337052662041377 0.316637311034378 0.12768360784343 -0.0293025419760333 0.011533333328731; - 0.00673504382217329 0.127802773462876 0.21400311568839 0.313312558886712 0.332320021608814 0.355738586947393 0.289302267356911 0.240342829317707 0.0668704675171058 -0.00673504382217329; - 0.0123456790123457 0.116567456572037 0.225284323338104 0.301940035273369 0.343862505804144 0.343862505804144 0.301940035273369 0.225284323338104 0.116567456572037 0.0123456790123457 - ] -const diffmat = [ - -27.1666666666667 33.1634374775264 -8.54863217041303 4 -2.42027662546121 1.70408819104185 -1.33333333333333 1.13247433143179 -1.03109120412576 0.5; - -8.29085936938159 4.01654328417507 5.75877048314363 -2.27431608520652 1.30540728933228 -0.898197570222574 0.694592710667722 -0.586256827714545 0.532088886237956 -0.257772801031441; - 2.13715804260326 -5.75877048314363 0.927019729872654 3.75877048314364 -1.68805925749197 1.06417777247591 -0.789861687269397 0.652703644666139 -0.586256827714545 0.283118582857949; - -1 2.27431608520652 -3.75877048314364 0.333333333333335 3.06417777247591 -1.48445439793712 1 -0.789861687269397 0.694592710667722 -0.333333333333333; - 0.605069156365302 -1.30540728933228 1.68805925749197 -3.06417777247591 0.0895235543024196 2.87938524157182 -1.48445439793712 1.06417777247591 -0.898197570222574 0.426022047760462; - -0.426022047760462 0.898197570222574 -1.06417777247591 1.48445439793712 -2.87938524157182 -0.0895235543024196 3.06417777247591 -1.68805925749197 1.30540728933228 -0.605069156365302; - 0.333333333333333 -0.694592710667722 0.789861687269397 -1 1.48445439793712 -3.06417777247591 -0.333333333333335 3.75877048314364 -2.27431608520652 1; - -0.283118582857949 0.586256827714545 -0.652703644666139 0.789861687269397 -1.06417777247591 1.68805925749197 -3.75877048314364 -0.927019729872654 5.75877048314363 -2.13715804260326; - 0.257772801031441 -0.532088886237956 0.586256827714545 -0.694592710667722 0.898197570222574 -1.30540728933228 2.27431608520652 -5.75877048314363 -4.01654328417507 8.29085936938159; - -0.5 1.03109120412576 -1.13247433143179 1.33333333333333 -1.70408819104185 2.42027662546121 -4 8.54863217041303 -33.1634374775264 27.1666666666667 - ] \ No newline at end of file +end \ No newline at end of file From f032ad081e728252763fa42003ca74d5567ba85b Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Tue, 5 Dec 2023 15:43:50 -0500 Subject: [PATCH 5/8] added test to cover barycentric interp, minor renaming of new functions --- src/gaussjacobi.jl | 6 +++--- test/test_gaussjacobi.jl | 12 ++++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index a7b6adb..cfa1ebf 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -377,7 +377,7 @@ function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) jk2 = .125*a8-(μ-1)./a8 - 4*(μ-1)*(7*μ-31)/3 ./ a8.^3 - 32*(μ-1)*(83*μ.^2-982*μ+3779)/15 ./ a8.^5 - 64*(μ-1)*(6949*μ^3-153855*μ^2+1585743*μ-6277237)/105 ./ a8.^7 - jk = [jk; jk2] + jk = append!(jk, jk2) end jk = real(jk[1:npts]) @@ -427,7 +427,7 @@ function feval_asy2(n::Integer, α::Float64, β::Float64, t::AbstractVector, hig Ja = besselj.(α, rho*t) Jb = besselj.(α + 1, rho*t) Jbb = besselj.(α + 1, rho2*t) - Jab = besselTaylor(-t, rho*t, α) + Jab = bessel_taylor(-t, rho*t, α) # Evaluate functions for recursive definition of coefficients: gt = A*(cot.(t/2) .- (2 ./ t)) .- B*tan.(t/2) @@ -628,7 +628,7 @@ function bary(x::AbstractVector, fvals::AbstractVector, xk::AbstractVector, vk:: return fx end -function besselTaylor(t::AbstractVector, z::AbstractVector, α::Float64) +function bessel_taylor(t::AbstractVector, z::AbstractVector, α::Float64) # Accurate evaluation of Bessel function J_A for asy2. (See [2].) # evaluates J_A(Z+T) by a Taylor series expansion about Z. diff --git a/test/test_gaussjacobi.jl b/test/test_gaussjacobi.jl index 6a5ce63..6ecfe7d 100644 --- a/test/test_gaussjacobi.jl +++ b/test/test_gaussjacobi.jl @@ -232,11 +232,23 @@ end @testset "boundary asymptotics were previously not implemented (#58,#130)" begin + # check example from issue #58 _, w = gaussjacobi(2^10, 0.25, 0) @test abs(w[end] - 3.607554904604311e-07) / 3.607554904604311e-07 < tol_w + # check example from issue #130 _, w = gaussjacobi(2^16, -0.9, 0) @test abs(w[end] - 1.223027243345722) < tol_w @test abs(sum(w) - 2^(-0.9+1)/(-0.9+1)) < tol + + # double check NaN fixing in barycentric interpolation implementation + x = range(0, stop=1, length=5) + N = 10 + fvals = ones(N) + xk = .5*(sin.(pi*(-(N-1):2:(N-1))/(2*(N-1))) .+ 1) + vk = [.5; ones(N-1)] + vk[2:2:end] .= -1 + vk[end] *= .5 + @test norm(FastGaussQuadrature.bary(x, fvals, xk, vk) - ones(5)) < tol end end \ No newline at end of file From 874a93c78253df486b54442837ca15e938d9ce03 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Wed, 6 Dec 2023 11:34:57 -0500 Subject: [PATCH 6/8] remove unreached lines for testing coverage --- src/gaussjacobi.jl | 13 +------------ test/test_gaussjacobi.jl | 6 +++--- 2 files changed, 4 insertions(+), 15 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index cfa1ebf..d93492d 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -368,18 +368,7 @@ function asy2(n::Integer, α::Float64, β::Float64, npts::Integer) # Algorithm for computing nodes and weights near the boundary. # Use Newton iterations to find the first few Bessel roots: - smallK = min(30, npts) - jk = approx_besselroots(α, smallK) - # Use asy formula for larger ones (See NIST 10.21.19, Olver 1974 p247) - if (npts > smallK) - μ = 4*α^2 - a8 = 8*(transpose(length(jk)+1:npts)+.5*α-.25)*pi - jk2 = .125*a8-(μ-1)./a8 - 4*(μ-1)*(7*μ-31)/3 ./ a8.^3 - - 32*(μ-1)*(83*μ.^2-982*μ+3779)/15 ./ a8.^5 - - 64*(μ-1)*(6949*μ^3-153855*μ^2+1585743*μ-6277237)/105 ./ a8.^7 - jk = append!(jk, jk2) - end - jk = real(jk[1:npts]) + jk = approx_besselroots(α, npts) # Approximate roots via asymptotic formula: (see Olver 1974, NIST, 18.16.8) phik = jk/(n + .5*(α + β + 1)) diff --git a/test/test_gaussjacobi.jl b/test/test_gaussjacobi.jl index 6ecfe7d..c777e21 100644 --- a/test/test_gaussjacobi.jl +++ b/test/test_gaussjacobi.jl @@ -234,11 +234,11 @@ @testset "boundary asymptotics were previously not implemented (#58,#130)" begin # check example from issue #58 _, w = gaussjacobi(2^10, 0.25, 0) - @test abs(w[end] - 3.607554904604311e-07) / 3.607554904604311e-07 < tol_w + @test abs(w[end] - 3.607554904604311e-07) < tol # check example from issue #130 _, w = gaussjacobi(2^16, -0.9, 0) - @test abs(w[end] - 1.223027243345722) < tol_w + @test abs(w[end] - 1.223027243345722) < tol @test abs(sum(w) - 2^(-0.9+1)/(-0.9+1)) < tol # double check NaN fixing in barycentric interpolation implementation @@ -249,6 +249,6 @@ vk = [.5; ones(N-1)] vk[2:2:end] .= -1 vk[end] *= .5 - @test norm(FastGaussQuadrature.bary(x, fvals, xk, vk) - ones(5)) < tol + @test FastGaussQuadrature.bary(x, fvals, xk, vk) ≈ ones(5) rtol=1e-16 end end \ No newline at end of file From 59802d553b4d0d24071184b9b348407c50383335 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Wed, 6 Dec 2023 20:21:24 -0500 Subject: [PATCH 7/8] remove innerjacobi_rec which has been replaced by feval_asy2 --- src/gaussjacobi.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/gaussjacobi.jl b/src/gaussjacobi.jl index d93492d..02a69d2 100644 --- a/src/gaussjacobi.jl +++ b/src/gaussjacobi.jl @@ -151,15 +151,6 @@ function innerjacobi_rec!(n, x, α::T, β::T, P, PP) where {T <: AbstractFloat} nothing end -function innerjacobi_rec(n, x, α::T, β::T) where {T <: AbstractFloat} - # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE. - N = length(x) - P = Array{T}(undef,N) - PP = Array{T}(undef,N) - innerjacobi_rec!(n, x, α, β, P, PP) - return P, PP -end - function weightsConstant(n, α, β) # Compute the constant for weights: M = min(20, n - 1) From 5a86bc13a5d814c40a06e6accd0b46b125837698 Mon Sep 17 00:00:00 2001 From: Paul Beckman Date: Thu, 7 Dec 2023 10:27:10 -0500 Subject: [PATCH 8/8] bump version to 1.0.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ab03eca..e7fef40 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastGaussQuadrature" uuid = "442a2c76-b920-505d-bb47-c5924d526838" -version = "1.0.0" +version = "1.0.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"