From f6d18213341545edc947935d84ad82624adea274 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 5 May 2025 20:41:21 +0100 Subject: [PATCH 1/3] Stieltjes on a square --- src/MultivariateSingularIntegrals.jl | 184 +-------------------------- src/logkernelsquare.jl | 181 ++++++++++++++++++++++++++ src/stieltjessquare.jl | 52 ++++++++ test/runtests.jl | 38 +++--- 4 files changed, 258 insertions(+), 197 deletions(-) create mode 100644 src/logkernelsquare.jl create mode 100644 src/stieltjessquare.jl diff --git a/src/MultivariateSingularIntegrals.jl b/src/MultivariateSingularIntegrals.jl index 1f5f088..ec21269 100644 --- a/src/MultivariateSingularIntegrals.jl +++ b/src/MultivariateSingularIntegrals.jl @@ -2,194 +2,16 @@ module MultivariateSingularIntegrals using LinearAlgebra, ClassicalOrthogonalPolynomials, SingularIntegrals, FillArrays import Base: size -export logkernelsquare +export logkernelsquare, stieltjessquare -""" - zlog(z) - -implements ``z*log(z)``. -""" -zlog(z) = iszero(z) ? zero(z) : z*log(z) - -""" - zlogm(z) - -implements ``z*log(z)`` but taking the other choice on the branch cut. -""" -zlogm(z) = iszero(imag(z)) && real(z) < 0 ? -zlog(abs(z)) - im*π*z : zlog(z) - -L0(z) = zlog(1 + complex(z)) - zlog(complex(z)-1) - 2one(z) -L1(z, r0=L0(z)) = (z+1)*r0/2 + 1 - zlog(complex(z)+1) - -function m_const(k, z) - (x,y) = reim(z) - T = float(real(typeof(z))) - im*convert(T,π) * if k == 0 - if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently - convert(T, 1) - elseif y ≤ -1 # x ≤ 0 - convert(T, -3) - else # x ≤ 0 and -1 < y < 1 - ((1 + y) - 3*(1-y))/2 - end - elseif -1 ≤ y ≤ 1 && x ≤ 0 # k ≠ 0 - -Weighted(Jacobi{T}(1,1))[y,k]/k - else - zero(T) - end -end - -M0(z) = L0(-im*float(z)) + m_const(0, float(z)) -M1(z, r0=L0(-im*z)) = L1(-im*float(z), r0) + m_const(1, float(z)) - - -""" - L₀₀(z) = ∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))dsdt -""" -L₀₀(z) = (1-z)*M0(z-1) + im*M1(z-1) + (1+z)*M0(z+1) - im*M1(z+1) - 4 - -function m_const_vec(n, z) - (x,y) = reim(z) - T = float(real(typeof(z))) - if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently - [im*convert(T,π); Zeros{T}(n-1)] - elseif y ≤ -1 - [-3im*convert(T,π); Zeros{T}(n-1)] - else # x ≤ 0 && -1 ≤ y ≤ 1 - r = Weighted(Jacobi{T}(1,1))[y,1:n-1] - [((1 + y) - 3*(1-y))/2*im*convert(T,π); -im .* convert(T,π) .* r ./ (1:(n-1))] - end -end function imlogkernel_vec(n, z) T = float(real(typeof(z))) transpose(complexlogkernel(Legendre{T}(), -im*float(z)))[1:n] + m_const_vec(n, float(z)) end - - -function rec_rhs_1!(F::AbstractMatrix{T}, z) where T - m,n = size(F) - x,y = reim(z) - πT = convert(T, π) - if x < -1 && -1 < y < 1 - C_y = Ultraspherical{T}(-1/2)[y,2:n+1] - F[1,:] .= (-4im*πT * x) .* C_y - elseif x < 1 && -1 < y < 1 - C_y = Ultraspherical{T}(-1/2)[y,2:n+1] - F[1,:] .= (-2im*πT * (x-1)) .* C_y - end - - F[1,1] += zlog(z-1-im) + zlogm(z-1+im) + zlog(z+1-im) + zlogm(z+1+im) - F[1,2] -= 4im/convert(T,3) - - M₋ = imlogkernel_vec(n+1, z-1) - M₊ = imlogkernel_vec(n+1, z+1) - - F[1,1] += im*(M₊[2] + M₋[2]) - for j = 1:n-1 - F[1,j+1] += im*(M₊[j+2] + M₋[j+2] - M₊[j] - M₋[j])/(2j+1) - end - F[2,1] -= 4/convert(T,3) - F -end - - - -function rec_rhs_2!(F::AbstractMatrix{T}, z) where T - m,n = size(F) - x,y = reim(z) - πT = convert(T, π) - if -1 < x < 1 && -1 ≤ y < 1 - C_x = Ultraspherical{T}(-3/2)[x,3:m+2] - C_y = Ultraspherical{T}(-1/2)[y,2:n+1] - F .= (2im*πT) .* (C_x .* C_y') ./ 3 - F[1,:] .-= (2im*πT) .* x .* C_y - F[2,:] .+= (2im*πT/3) .* C_y - elseif x ≤ -1 && -1 ≤ y < 1 - fill!(F, zero(T)) - C_y = Ultraspherical{T}(-1/2)[y,2:n+1] - F[1,:] .= (-4im*πT) .* x .* C_y - F[2,:] .= (4im*πT/3) .* C_y - else - fill!(F, zero(T)) - end - - L₋ = complexlogkernel(Legendre{T}(), z-im)[1:m+1] - L₊ = complexlogkernel(Legendre{T}(), z+im)[1:m+1] - - F[1,1] += L₋[2] + L₊[2] + zlog(z-im-1) + zlog(z-im+1) + zlog(z+im-1) + zlog(z+im+1) - F[1,2] -= 4im/convert(T,3) - for k = 1:m-1 - F[k+1,1] += (L₊[k+2] + L₋[k+2] - L₊[k] - L₋[k])/(2k+1) - end - F[2,1] -= 4/convert(T,3) - F -end - - -function logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2) - A[2,1] = z * A[1,1]/3 + (F_2[1,1] - 2F_1[1,1])/3 - for k = 1:size(A,1)-2 - A[k+2,1] = (2k+1) * z * A[k+1,1]/(k+3) - (k-2)*A[k,1]/(k+3) + (2k+1) * (F_2[k+1,1] - 2F_1[k+1,1])/(k+3) - end - A -end - -function logkernelsquare_populatefirstrow!(A, z, F_1, F_2) - A[1,2] = z*A[1,1]/(3im) + (F_1[1,1]-2F_2[1,1])/(3im) - for j = 1:size(A,2)-2 - A[1,j+2] = -im * (2j+1) * z * A[1,j+1]/(j+3) - (j-2)*A[1,j]/(j+3) - im * (2j+1) * (F_1[1,j+1] - 2F_2[1,j+1])/(j+3) - end - A -end - - -""" - logkernelsquare(z, n) - -computes the matrix with entries ``∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))P_k(s)P_j(t)dsdt`` up to total degree ``n``. -The bottom right of the returned matrix is zero. For a square truncation compute `logkernelsquare(z,2n-1)[1:n,1:n]`. -""" - -function logkernelsquare(z, n) - T = complex(float(eltype(z))) - logkernelsquare!(zeros(T,n,n), z, zeros(T,n,n), zeros(T,n,n)) -end - - -function logkernelsquare!(A::AbstractMatrix{T}, z, F_1, F_2) where T - m,n = size(A) - @assert m == n - rec_rhs_1!(F_1, z) - rec_rhs_2!(F_2, z) - A[1,1] = L₀₀(z) - logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2) - logkernelsquare_populatefirstrow!(A, z, F_1, F_2) - - # F = F_1 # reuse the memory - F = F_2 .- F_1 - - # 2nd row/column - - for k = 1:m-2 - A[k+1,2] = im*(F[k+1,1] + (A[k,1] - A[k+2,1])/(2k+1)) - end - - for j = 2:n-2 - A[2,j+1] = F[1,j+1] + im*(A[1,j+2] - A[1,j])/(2j+1) - end - - for ℓ = 1:((n-1)÷2-1) - for k = ℓ+1:n-(ℓ+2) - A[k+1,ℓ+2] = (2ℓ+1)*im*(F[k+1,ℓ+1] + (A[k,ℓ+1] - A[k+2,ℓ+1])/(2k+1)) + A[k+1,ℓ] - end - for j = ℓ+2:n-(ℓ+2) - A[ℓ+2,j+1] = (2ℓ+1) * (F[ℓ+1,j+1] + im*(A[ℓ+1,j+2] - A[ℓ+1,j])/(2j+1)) + A[ℓ,j+1] - end - end - A -end +include("stieltjessquare.jl") +include("logkernelsquare.jl") end # module MultivariateSingularIntegrals diff --git a/src/logkernelsquare.jl b/src/logkernelsquare.jl new file mode 100644 index 0000000..1a0e15b --- /dev/null +++ b/src/logkernelsquare.jl @@ -0,0 +1,181 @@ +""" + zlog(z) + +implements ``z*log(z)``. +""" +zlog(z) = iszero(z) ? zero(z) : z*log(z) + +""" + zlogm(z) + +implements ``z*log(z)`` but taking the other choice on the branch cut. +""" +zlogm(z) = iszero(imag(z)) && real(z) < 0 ? -zlog(abs(z)) - im*π*z : zlog(z) + +L0(z) = zlog(1 + complex(z)) - zlog(complex(z)-1) - 2one(z) +L1(z, r0=L0(z)) = (z+1)*r0/2 + 1 - zlog(complex(z)+1) + +function m_const(k, z) + (x,y) = reim(z) + T = float(real(typeof(z))) + im*convert(T,π) * if k == 0 + if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently + convert(T, 1) + elseif y ≤ -1 # x ≤ 0 + convert(T, -3) + else # x ≤ 0 and -1 < y < 1 + ((1 + y) - 3*(1-y))/2 + end + elseif -1 ≤ y ≤ 1 && x ≤ 0 # k ≠ 0 + -Weighted(Jacobi{T}(1,1))[y,k]/k + else + zero(T) + end +end + +M0(z) = L0(-im*float(z)) + m_const(0, float(z)) +M1(z, r0=L0(-im*z)) = L1(-im*float(z), r0) + m_const(1, float(z)) + + +""" + L₀₀(z) = ∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))dsdt +""" +L₀₀(z) = (1-z)*M0(z-1) + im*M1(z-1) + (1+z)*M0(z+1) - im*M1(z+1) - 4 + +function m_const_vec(n, z) + (x,y) = reim(z) + T = float(real(typeof(z))) + if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently + [im*convert(T,π); Zeros{T}(n-1)] + elseif y ≤ -1 + [-3im*convert(T,π); Zeros{T}(n-1)] + else # x ≤ 0 && -1 ≤ y ≤ 1 + r = Weighted(Jacobi{T}(1,1))[y,1:n-1] + [((1 + y) - 3*(1-y))/2*im*convert(T,π); -im .* convert(T,π) .* r ./ (1:(n-1))] + end +end + + + +function rec_rhs_1!(F::AbstractMatrix{T}, z) where T + m,n = size(F) + x,y = reim(z) + πT = convert(T, π) + if x < -1 && -1 < y < 1 + C_y = Ultraspherical{T}(-1/2)[y,2:n+1] + F[1,:] .= (-4im*πT * x) .* C_y + elseif x < 1 && -1 < y < 1 + C_y = Ultraspherical{T}(-1/2)[y,2:n+1] + F[1,:] .= (-2im*πT * (x-1)) .* C_y + end + + F[1,1] += zlog(z-1-im) + zlogm(z-1+im) + zlog(z+1-im) + zlogm(z+1+im) + F[1,2] -= 4im/convert(T,3) + + M₋ = imlogkernel_vec(n+1, z-1) + M₊ = imlogkernel_vec(n+1, z+1) + + F[1,1] += im*(M₊[2] + M₋[2]) + for j = 1:n-1 + F[1,j+1] += im*(M₊[j+2] + M₋[j+2] - M₊[j] - M₋[j])/(2j+1) + end + F[2,1] -= 4/convert(T,3) + F +end + + + +function rec_rhs_2!(F::AbstractMatrix{T}, z) where T + m,n = size(F) + x,y = reim(z) + πT = convert(T, π) + if -1 < x < 1 && -1 ≤ y < 1 + C_x = Ultraspherical{T}(-3/2)[x,3:m+2] + C_y = Ultraspherical{T}(-1/2)[y,2:n+1] + F .= (2im*πT) .* (C_x .* C_y') ./ 3 + F[1,:] .-= (2im*πT) .* x .* C_y + F[2,:] .+= (2im*πT/3) .* C_y + elseif x ≤ -1 && -1 ≤ y < 1 + fill!(F, zero(T)) + C_y = Ultraspherical{T}(-1/2)[y,2:n+1] + F[1,:] .= (-4im*πT) .* x .* C_y + F[2,:] .= (4im*πT/3) .* C_y + else + fill!(F, zero(T)) + end + + L₋ = complexlogkernel(Legendre{T}(), z-im)[1:m+1] + L₊ = complexlogkernel(Legendre{T}(), z+im)[1:m+1] + + F[1,1] += L₋[2] + L₊[2] + zlog(z-im-1) + zlog(z-im+1) + zlog(z+im-1) + zlog(z+im+1) + F[1,2] -= 4im/convert(T,3) + for k = 1:m-1 + F[k+1,1] += (L₊[k+2] + L₋[k+2] - L₊[k] - L₋[k])/(2k+1) + end + F[2,1] -= 4/convert(T,3) + F +end + + +function logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2) + A[2,1] = z * A[1,1]/3 + (F_2[1,1] - 2F_1[1,1])/3 + for k = 1:size(A,1)-2 + A[k+2,1] = (2k+1) * z * A[k+1,1]/(k+3) - (k-2)*A[k,1]/(k+3) + (2k+1) * (F_2[k+1,1] - 2F_1[k+1,1])/(k+3) + end + A +end + +function logkernelsquare_populatefirstrow!(A, z, F_1, F_2) + A[1,2] = z*A[1,1]/(3im) + (F_1[1,1]-2F_2[1,1])/(3im) + for j = 1:size(A,2)-2 + A[1,j+2] = -im * (2j+1) * z * A[1,j+1]/(j+3) - (j-2)*A[1,j]/(j+3) - im * (2j+1) * (F_1[1,j+1] - 2F_2[1,j+1])/(j+3) + end + A +end + + +""" + logkernelsquare(z, n) + +computes the matrix with entries ``∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))P_k(s)P_j(t)dsdt`` up to total degree ``n``. +The bottom right of the returned matrix is zero. For a square truncation compute `logkernelsquare(z,2n-1)[1:n,1:n]`. +""" + +function logkernelsquare(z, n) + T = complex(float(eltype(z))) + logkernelsquare!(zeros(T,n,n), z, zeros(T,n,n), zeros(T,n,n)) +end + + +function logkernelsquare!(A::AbstractMatrix{T}, z, F_1, F_2) where T + m,n = size(A) + @assert m == n + rec_rhs_1!(F_1, z) + rec_rhs_2!(F_2, z) + A[1,1] = L₀₀(z) + logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2) + logkernelsquare_populatefirstrow!(A, z, F_1, F_2) + + # F = F_1 # reuse the memory + F = F_2 .- F_1 + + # 2nd row/column + + for k = 1:m-2 + A[k+1,2] = im*(F[k+1,1] + (A[k,1] - A[k+2,1])/(2k+1)) + end + + for j = 2:n-2 + A[2,j+1] = F[1,j+1] + im*(A[1,j+2] - A[1,j])/(2j+1) + end + + for ℓ = 1:((n-1)÷2-1) + for k = ℓ+1:n-(ℓ+2) + A[k+1,ℓ+2] = (2ℓ+1)*im*(F[k+1,ℓ+1] + (A[k,ℓ+1] - A[k+2,ℓ+1])/(2k+1)) + A[k+1,ℓ] + end + for j = ℓ+2:n-(ℓ+2) + A[ℓ+2,j+1] = (2ℓ+1) * (F[ℓ+1,j+1] + im*(A[ℓ+1,j+2] - A[ℓ+1,j])/(2j+1)) + A[ℓ,j+1] + end + end + A +end \ No newline at end of file diff --git a/src/stieltjessquare.jl b/src/stieltjessquare.jl new file mode 100644 index 0000000..f41a84d --- /dev/null +++ b/src/stieltjessquare.jl @@ -0,0 +1,52 @@ +""" + stieltjessquare(z, n) + +computes the matrix with entries ``∫_{-1}^1∫_{-1}^1 P_k(s)P_j(t)/(z-(s+im*t)) dsdt`` up to total degree ``n``. +The bottom right of the returned matrix is zero. For a square truncation compute `stieltjessquare(z,2n-1)[1:n,1:n]`. +""" + +function stieltjessquare(z, n) + T = complex(float(eltype(z))) + stieltjessquare!(zeros(T,n,n), z) +end + +function stieltjessquare_populatefirstcolumn!(A, z) + M_p = imlogkernel_vec(n, z+1) + M_m = imlogkernel_vec(n, z-1) + A[:,1] .= M_p .- M_m +end + +function stieltjessquare_populatefirstrow!(A, z) + M_p = imlogkernel_vec(n, 1-im*z) + M_m = imlogkernel_vec(n, -1-im*z) + m = size(A,1) + A[1,:] .= (-1) .^ (1:m) .* im .* (M_p .- M_m) +end + + +function stieltjessquare!(A::AbstractMatrix{T}, z) where T + m,n = size(A) + @assert m == n + + stieltjessquare_populatefirstcolumn!(A, z, F_1, F_2) + stieltjessquare_populatefirstrow!(A, z, F_1, F_2) + + # 2nd row/column + for k = 1:m-2 + A[k+1,2] = im* (k*A[k,1]/(2k+1) + (k+1)*A[k+2,1]/(2k+1) - z*A[k+1,1]) + end + + for j = 2:n-2 + A[2,j+1] = z*A[1,j+1] - im*(j*A[1,j]/(2j+1) + (j+1)*A[1,k+2]/(2j+1)) + end + # remaining + for ℓ = 1:((n-1)÷2-1) + for k = ℓ+1:n-(ℓ+2) + A[k+1,ℓ+2] = (im*(2ℓ+1)*(k*A[k,ℓ+1]/(2k+1) + (k+1)*A[k+2,ℓ+1]/(2k+1) - z*A[k+1,ℓ+1]) - j*A[k+1,ℓ+1])/(ℓ+1) + end + for j = ℓ+2:n-(ℓ+2) + A[ℓ+2,j+1] = (z*(2k+1)*A[ℓ+1,j+1] - im*(2k+1)*(j*A[ℓ+1,j]/(2j+1) + (j+1)*A[ℓ+1,k+2]/(2j+1)) - k*A[k,ℓ+1])/(k+1) + end + end + A +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 1da5bd7..9d0b157 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,24 +1,30 @@ using MultivariateSingularIntegrals, QuadGK, ClassicalOrthogonalPolynomials, Test -@testset "accurate for low order" begin - Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im) - n = 6 - for z in Z - L = logkernelsquare(z,n) - for j = 0:n-1, k=0:n-j-1 - q = quadgk(s -> quadgk(t -> iszero(s+im*t) ? 0.0+0.0im : log(z-(s+im*t)) * legendrep(k,s) * legendrep(j,t), -1, 1, atol=1E-12)[1], -1, 1, atol=1E-12)[1] - @test q ≈ L[k+1,j+1] atol=1E-10 +Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im) +@testset "stieltjes" begin + +end + +@testset "logkernel" begin + @testset "accurate for low order" begin + n = 6 + for z in Z + L = logkernelsquare(z,n) + for j = 0:n-1, k=0:n-j-1 + q = quadgk(s -> quadgk(t -> iszero(s+im*t) ? 0.0+0.0im : log(z-(s+im*t)) * legendrep(k,s) * legendrep(j,t), -1, 1, atol=1E-12)[1], -1, 1, atol=1E-12)[1] + @test q ≈ L[k+1,j+1] atol=1E-10 + end end end -end -@testset "BigFloat" begin - setprecision(2000) do - for z in (big(2.0), big(2.0im), big(2.0+2im), big(2.0-2im), big(-2.0 - 2im)) - P = Legendre{BigFloat}() - n = 100 - M = Diagonal((P'P)[1:n,1:n]) - @test P[big(0.), 1:n]' * inv(M) * logkernelsquare(z, n) * inv(M) * P[big(0.), 1:n] ≈ log(z) atol=1E-40 + @testset "BigFloat" begin + setprecision(2000) do + for z in (big(2.0), big(2.0im), big(2.0+2im), big(2.0-2im), big(-2.0 - 2im)) + P = Legendre{BigFloat}() + n = 100 + M = Diagonal((P'P)[1:n,1:n]) + @test P[big(0.), 1:n]' * inv(M) * logkernelsquare(z, n) * inv(M) * P[big(0.), 1:n] ≈ log(z) atol=1E-40 + end end end end From 73c672330dfccebb055e07e0949bcbb2c6b83189 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 6 May 2025 10:24:57 +0100 Subject: [PATCH 2/3] Update runtests.jl --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9d0b157..04b3634 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using MultivariateSingularIntegrals, QuadGK, ClassicalOrthogonalPolynomials, Test +using MultivariateSingularIntegrals, QuadGK, ClassicalOrthogonalPolynomials, LinearAlgebra, Test Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im) @testset "stieltjes" begin From 3427a72aeddec4988dcd3727897b3249584c57a4 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 6 May 2025 11:24:11 +0100 Subject: [PATCH 3/3] Add tests and fix Stieltjes --- src/stieltjessquare.jl | 25 +++++++++++++------------ test/runtests.jl | 24 +++++++++++++++++++++++- 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/src/stieltjessquare.jl b/src/stieltjessquare.jl index f41a84d..236dc50 100644 --- a/src/stieltjessquare.jl +++ b/src/stieltjessquare.jl @@ -10,17 +10,18 @@ function stieltjessquare(z, n) stieltjessquare!(zeros(T,n,n), z) end -function stieltjessquare_populatefirstcolumn!(A, z) +function stieltjessquare_populatefirstrow!(A, z) + n = size(A,2) M_p = imlogkernel_vec(n, z+1) M_m = imlogkernel_vec(n, z-1) - A[:,1] .= M_p .- M_m + A[1,:] .= M_p .- M_m end -function stieltjessquare_populatefirstrow!(A, z) - M_p = imlogkernel_vec(n, 1-im*z) - M_m = imlogkernel_vec(n, -1-im*z) - m = size(A,1) - A[1,:] .= (-1) .^ (1:m) .* im .* (M_p .- M_m) +function stieltjessquare_populatefirstcolumn!(A, z) + m = size(A,2) + M_p = imlogkernel_vec(m, 1-im*z) + M_m = imlogkernel_vec(m, -1-im*z) + A[:,1] .= (-1) .^ (1:m) .* im .* (M_p .- M_m) end @@ -28,8 +29,8 @@ function stieltjessquare!(A::AbstractMatrix{T}, z) where T m,n = size(A) @assert m == n - stieltjessquare_populatefirstcolumn!(A, z, F_1, F_2) - stieltjessquare_populatefirstrow!(A, z, F_1, F_2) + stieltjessquare_populatefirstcolumn!(A, z) + stieltjessquare_populatefirstrow!(A, z) # 2nd row/column for k = 1:m-2 @@ -37,15 +38,15 @@ function stieltjessquare!(A::AbstractMatrix{T}, z) where T end for j = 2:n-2 - A[2,j+1] = z*A[1,j+1] - im*(j*A[1,j]/(2j+1) + (j+1)*A[1,k+2]/(2j+1)) + A[2,j+1] = z*A[1,j+1] - im*(j*A[1,j]/(2j+1) + (j+1)*A[1,j+2]/(2j+1)) end # remaining for ℓ = 1:((n-1)÷2-1) for k = ℓ+1:n-(ℓ+2) - A[k+1,ℓ+2] = (im*(2ℓ+1)*(k*A[k,ℓ+1]/(2k+1) + (k+1)*A[k+2,ℓ+1]/(2k+1) - z*A[k+1,ℓ+1]) - j*A[k+1,ℓ+1])/(ℓ+1) + A[k+1,ℓ+2] = (-im * (2ℓ+1)* z * A[k+1,ℓ+1] + im*(2ℓ+1)*(k*A[k,ℓ+1] + (k+1)*A[k+2,ℓ+1])/(2k+1) - ℓ*A[k+1,ℓ])/(ℓ+1) end for j = ℓ+2:n-(ℓ+2) - A[ℓ+2,j+1] = (z*(2k+1)*A[ℓ+1,j+1] - im*(2k+1)*(j*A[ℓ+1,j]/(2j+1) + (j+1)*A[ℓ+1,k+2]/(2j+1)) - k*A[k,ℓ+1])/(k+1) + A[ℓ+2,j+1] = (z*(2ℓ+1)*A[ℓ+1,j+1] - im*(2ℓ+1)*(j*A[ℓ+1,j] + (j+1)*A[ℓ+1,j+2])/(2j+1) - ℓ*A[ℓ,j+1])/(ℓ+1) end end A diff --git a/test/runtests.jl b/test/runtests.jl index 04b3634..cc24274 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,30 @@ using MultivariateSingularIntegrals, QuadGK, ClassicalOrthogonalPolynomials, LinearAlgebra, Test + +Z̃ = (2.0, 2.0im, 2.0+2im, 2.0-2im, -2.0 - 2im) Z = (2, 2im, 2+2im, 2-2im, -2 - 2im, -1.001-1.5im, -0.999-1.5im, -0.5-2im, -1-1.5im, -2+2im, -2, 0.1+0.2im, 0.1+im, 0.1-im, 1+0.1im, -1+0.1im, 1+im, -1+im, -1-im, 1-im) @testset "stieltjes" begin - + @testset "accurate for low order" begin + n = 6 + for z in Z̃ + L = stieltjessquare(z,n) + for j = 0:n-1, k=0:n-j-1 + q = quadgk(s -> quadgk(t -> iszero(s+im*t) ? 0.0+0.0im : inv(z-(s+im*t)) * legendrep(k,s) * legendrep(j,t), -1, 1, atol=1E-12)[1], -1, 1, atol=1E-12)[1] + @test q ≈ L[k+1,j+1] atol=1E-10 + end + end + end + + @testset "BigFloat" begin + setprecision(2000) do + for z in big.(Z̃) + P = Legendre{BigFloat}() + n = 100 + M = Diagonal((P'P)[1:n,1:n]) + @test P[big(0.), 1:n]' * inv(M) * stieltjessquare(z, n) * inv(M) * P[big(0.), 1:n] ≈ inv(z) atol=1E-40 + end + end + end end @testset "logkernel" begin