Skip to content

Commit d00728a

Browse files
dlfivefiftydlfivefifty
authored
Stieltjes on a square (#1)
* Stieltjes on a square * Update runtests.jl * Add tests and fix Stieltjes
1 parent 4e23613 commit d00728a

File tree

4 files changed

+281
-197
lines changed

4 files changed

+281
-197
lines changed

‎src/MultivariateSingularIntegrals.jl

Lines changed: 3 additions & 181 deletions
Original file line numberDiff line numberDiff line change
@@ -2,194 +2,16 @@ module MultivariateSingularIntegrals
22
using LinearAlgebra, ClassicalOrthogonalPolynomials, SingularIntegrals, FillArrays
33
import Base: size
44

5-
export logkernelsquare
5+
export logkernelsquare, stieltjessquare
66

7-
"""
8-
zlog(z)
9-
10-
implements ``z*log(z)``.
11-
"""
12-
zlog(z) = iszero(z) ? zero(z) : z*log(z)
13-
14-
"""
15-
zlogm(z)
16-
17-
implements ``z*log(z)`` but taking the other choice on the branch cut.
18-
"""
19-
zlogm(z) = iszero(imag(z)) && real(z) < 0 ? -zlog(abs(z)) - im*π*z : zlog(z)
20-
21-
L0(z) = zlog(1 + complex(z)) - zlog(complex(z)-1) - 2one(z)
22-
L1(z, r0=L0(z)) = (z+1)*r0/2 + 1 - zlog(complex(z)+1)
23-
24-
function m_const(k, z)
25-
(x,y) = reim(z)
26-
T = float(real(typeof(z)))
27-
im*convert(T,π) * if k == 0
28-
if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently
29-
convert(T, 1)
30-
elseif y ≤ -1 # x ≤ 0
31-
convert(T, -3)
32-
else # x ≤ 0 and -1 < y < 1
33-
((1 + y) - 3*(1-y))/2
34-
end
35-
elseif -1 ≤ y ≤ 1 && x ≤ 0 # k ≠ 0
36-
-Weighted(Jacobi{T}(1,1))[y,k]/k
37-
else
38-
zero(T)
39-
end
40-
end
41-
42-
M0(z) = L0(-im*float(z)) + m_const(0, float(z))
43-
M1(z, r0=L0(-im*z)) = L1(-im*float(z), r0) + m_const(1, float(z))
44-
45-
46-
"""
47-
L₀₀(z) = ∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))dsdt
48-
"""
49-
L₀₀(z) = (1-z)*M0(z-1) + im*M1(z-1) + (1+z)*M0(z+1) - im*M1(z+1) - 4
50-
51-
function m_const_vec(n, z)
52-
(x,y) = reim(z)
53-
T = float(real(typeof(z)))
54-
if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently
55-
[im*convert(T,π); Zeros{T}(n-1)]
56-
elseif y ≤ -1
57-
[-3im*convert(T,π); Zeros{T}(n-1)]
58-
else # x ≤ 0 && -1 ≤ y ≤ 1
59-
r = Weighted(Jacobi{T}(1,1))[y,1:n-1]
60-
[((1 + y) - 3*(1-y))/2*im*convert(T,π); -im .* convert(T,π) .* r ./ (1:(n-1))]
61-
end
62-
end
637

648
function imlogkernel_vec(n, z)
659
T = float(real(typeof(z)))
6610
transpose(complexlogkernel(Legendre{T}(), -im*float(z)))[1:n] + m_const_vec(n, float(z))
6711
end
6812

69-
70-
71-
function rec_rhs_1!(F::AbstractMatrix{T}, z) where T
72-
m,n = size(F)
73-
x,y = reim(z)
74-
πT = convert(T, π)
75-
if x < -1 && -1 < y < 1
76-
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
77-
F[1,:] .= (-4im*πT * x) .* C_y
78-
elseif x < 1 && -1 < y < 1
79-
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
80-
F[1,:] .= (-2im*πT * (x-1)) .* C_y
81-
end
82-
83-
F[1,1] += zlog(z-1-im) + zlogm(z-1+im) + zlog(z+1-im) + zlogm(z+1+im)
84-
F[1,2] -= 4im/convert(T,3)
85-
86-
M₋ = imlogkernel_vec(n+1, z-1)
87-
M₊ = imlogkernel_vec(n+1, z+1)
88-
89-
F[1,1] += im*(M₊[2] + M₋[2])
90-
for j = 1:n-1
91-
F[1,j+1] += im*(M₊[j+2] + M₋[j+2] - M₊[j] - M₋[j])/(2j+1)
92-
end
93-
F[2,1] -= 4/convert(T,3)
94-
F
95-
end
96-
97-
98-
99-
function rec_rhs_2!(F::AbstractMatrix{T}, z) where T
100-
m,n = size(F)
101-
x,y = reim(z)
102-
πT = convert(T, π)
103-
if -1 < x < 1 && -1 ≤ y < 1
104-
C_x = Ultraspherical{T}(-3/2)[x,3:m+2]
105-
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
106-
F .= (2im*πT) .* (C_x .* C_y') ./ 3
107-
F[1,:] .-= (2im*πT) .* x .* C_y
108-
F[2,:] .+= (2im*πT/3) .* C_y
109-
elseif x ≤ -1 && -1 ≤ y < 1
110-
fill!(F, zero(T))
111-
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
112-
F[1,:] .= (-4im*πT) .* x .* C_y
113-
F[2,:] .= (4im*πT/3) .* C_y
114-
else
115-
fill!(F, zero(T))
116-
end
117-
118-
L₋ = complexlogkernel(Legendre{T}(), z-im)[1:m+1]
119-
L₊ = complexlogkernel(Legendre{T}(), z+im)[1:m+1]
120-
121-
F[1,1] += L₋[2] + L₊[2] + zlog(z-im-1) + zlog(z-im+1) + zlog(z+im-1) + zlog(z+im+1)
122-
F[1,2] -= 4im/convert(T,3)
123-
for k = 1:m-1
124-
F[k+1,1] += (L₊[k+2] + L₋[k+2] - L₊[k] - L₋[k])/(2k+1)
125-
end
126-
F[2,1] -= 4/convert(T,3)
127-
F
128-
end
129-
130-
131-
function logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2)
132-
A[2,1] = z * A[1,1]/3 + (F_2[1,1] - 2F_1[1,1])/3
133-
for k = 1:size(A,1)-2
134-
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)
135-
end
136-
A
137-
end
138-
139-
function logkernelsquare_populatefirstrow!(A, z, F_1, F_2)
140-
A[1,2] = z*A[1,1]/(3im) + (F_1[1,1]-2F_2[1,1])/(3im)
141-
for j = 1:size(A,2)-2
142-
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)
143-
end
144-
A
145-
end
146-
147-
148-
"""
149-
logkernelsquare(z, n)
150-
151-
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``.
152-
The bottom right of the returned matrix is zero. For a square truncation compute `logkernelsquare(z,2n-1)[1:n,1:n]`.
153-
"""
154-
155-
function logkernelsquare(z, n)
156-
T = complex(float(eltype(z)))
157-
logkernelsquare!(zeros(T,n,n), z, zeros(T,n,n), zeros(T,n,n))
158-
end
159-
160-
161-
function logkernelsquare!(A::AbstractMatrix{T}, z, F_1, F_2) where T
162-
m,n = size(A)
163-
@assert m == n
164-
rec_rhs_1!(F_1, z)
165-
rec_rhs_2!(F_2, z)
166-
A[1,1] = L₀₀(z)
167-
logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2)
168-
logkernelsquare_populatefirstrow!(A, z, F_1, F_2)
169-
170-
# F = F_1 # reuse the memory
171-
F = F_2 .- F_1
172-
173-
# 2nd row/column
174-
175-
for k = 1:m-2
176-
A[k+1,2] = im*(F[k+1,1] + (A[k,1] - A[k+2,1])/(2k+1))
177-
end
178-
179-
for j = 2:n-2
180-
A[2,j+1] = F[1,j+1] + im*(A[1,j+2] - A[1,j])/(2j+1)
181-
end
182-
183-
for ℓ = 1:((n-1)÷2-1)
184-
for k = ℓ+1:n-(ℓ+2)
185-
A[k+1,ℓ+2] = (2ℓ+1)*im*(F[k+1,ℓ+1] + (A[k,ℓ+1] - A[k+2,ℓ+1])/(2k+1)) + A[k+1,ℓ]
186-
end
187-
for j = ℓ+2:n-(ℓ+2)
188-
A[ℓ+2,j+1] = (2ℓ+1) * (F[ℓ+1,j+1] + im*(A[ℓ+1,j+2] - A[ℓ+1,j])/(2j+1)) + A[ℓ,j+1]
189-
end
190-
end
191-
A
192-
end
13+
include("stieltjessquare.jl")
14+
include("logkernelsquare.jl")
19315

19416

19517
end # module MultivariateSingularIntegrals

‎src/logkernelsquare.jl

Lines changed: 181 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,181 @@
1+
"""
2+
zlog(z)
3+
4+
implements ``z*log(z)``.
5+
"""
6+
zlog(z) = iszero(z) ? zero(z) : z*log(z)
7+
8+
"""
9+
zlogm(z)
10+
11+
implements ``z*log(z)`` but taking the other choice on the branch cut.
12+
"""
13+
zlogm(z) = iszero(imag(z)) && real(z) < 0 ? -zlog(abs(z)) - im*π*z : zlog(z)
14+
15+
L0(z) = zlog(1 + complex(z)) - zlog(complex(z)-1) - 2one(z)
16+
L1(z, r0=L0(z)) = (z+1)*r0/2 + 1 - zlog(complex(z)+1)
17+
18+
function m_const(k, z)
19+
(x,y) = reim(z)
20+
T = float(real(typeof(z)))
21+
im*convert(T,π) * if k == 0
22+
if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently
23+
convert(T, 1)
24+
elseif y ≤ -1 # x ≤ 0
25+
convert(T, -3)
26+
else # x ≤ 0 and -1 < y < 1
27+
((1 + y) - 3*(1-y))/2
28+
end
29+
elseif -1 ≤ y ≤ 1 && x ≤ 0 # k ≠ 0
30+
-Weighted(Jacobi{T}(1,1))[y,k]/k
31+
else
32+
zero(T)
33+
end
34+
end
35+
36+
M0(z) = L0(-im*float(z)) + m_const(0, float(z))
37+
M1(z, r0=L0(-im*z)) = L1(-im*float(z), r0) + m_const(1, float(z))
38+
39+
40+
"""
41+
L₀₀(z) = ∫_{-1}^1∫_{-1}^1 log(z-(s+im*t))dsdt
42+
"""
43+
L₀₀(z) = (1-z)*M0(z-1) + im*M1(z-1) + (1+z)*M0(z+1) - im*M1(z+1) - 4
44+
45+
function m_const_vec(n, z)
46+
(x,y) = reim(z)
47+
T = float(real(typeof(z)))
48+
if x > 0 || y ≥ 1 || (x == 0 && !signbit(x) && y < 0) # need to treat im*y+0.0 differently
49+
[im*convert(T,π); Zeros{T}(n-1)]
50+
elseif y ≤ -1
51+
[-3im*convert(T,π); Zeros{T}(n-1)]
52+
else # x ≤ 0 && -1 ≤ y ≤ 1
53+
r = Weighted(Jacobi{T}(1,1))[y,1:n-1]
54+
[((1 + y) - 3*(1-y))/2*im*convert(T,π); -im .* convert(T,π) .* r ./ (1:(n-1))]
55+
end
56+
end
57+
58+
59+
60+
function rec_rhs_1!(F::AbstractMatrix{T}, z) where T
61+
m,n = size(F)
62+
x,y = reim(z)
63+
πT = convert(T, π)
64+
if x < -1 && -1 < y < 1
65+
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
66+
F[1,:] .= (-4im*πT * x) .* C_y
67+
elseif x < 1 && -1 < y < 1
68+
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
69+
F[1,:] .= (-2im*πT * (x-1)) .* C_y
70+
end
71+
72+
F[1,1] += zlog(z-1-im) + zlogm(z-1+im) + zlog(z+1-im) + zlogm(z+1+im)
73+
F[1,2] -= 4im/convert(T,3)
74+
75+
M₋ = imlogkernel_vec(n+1, z-1)
76+
M₊ = imlogkernel_vec(n+1, z+1)
77+
78+
F[1,1] += im*(M₊[2] + M₋[2])
79+
for j = 1:n-1
80+
F[1,j+1] += im*(M₊[j+2] + M₋[j+2] - M₊[j] - M₋[j])/(2j+1)
81+
end
82+
F[2,1] -= 4/convert(T,3)
83+
F
84+
end
85+
86+
87+
88+
function rec_rhs_2!(F::AbstractMatrix{T}, z) where T
89+
m,n = size(F)
90+
x,y = reim(z)
91+
πT = convert(T, π)
92+
if -1 < x < 1 && -1 ≤ y < 1
93+
C_x = Ultraspherical{T}(-3/2)[x,3:m+2]
94+
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
95+
F .= (2im*πT) .* (C_x .* C_y') ./ 3
96+
F[1,:] .-= (2im*πT) .* x .* C_y
97+
F[2,:] .+= (2im*πT/3) .* C_y
98+
elseif x ≤ -1 && -1 ≤ y < 1
99+
fill!(F, zero(T))
100+
C_y = Ultraspherical{T}(-1/2)[y,2:n+1]
101+
F[1,:] .= (-4im*πT) .* x .* C_y
102+
F[2,:] .= (4im*πT/3) .* C_y
103+
else
104+
fill!(F, zero(T))
105+
end
106+
107+
L₋ = complexlogkernel(Legendre{T}(), z-im)[1:m+1]
108+
L₊ = complexlogkernel(Legendre{T}(), z+im)[1:m+1]
109+
110+
F[1,1] += L₋[2] + L₊[2] + zlog(z-im-1) + zlog(z-im+1) + zlog(z+im-1) + zlog(z+im+1)
111+
F[1,2] -= 4im/convert(T,3)
112+
for k = 1:m-1
113+
F[k+1,1] += (L₊[k+2] + L₋[k+2] - L₊[k] - L₋[k])/(2k+1)
114+
end
115+
F[2,1] -= 4/convert(T,3)
116+
F
117+
end
118+
119+
120+
function logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2)
121+
A[2,1] = z * A[1,1]/3 + (F_2[1,1] - 2F_1[1,1])/3
122+
for k = 1:size(A,1)-2
123+
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)
124+
end
125+
A
126+
end
127+
128+
function logkernelsquare_populatefirstrow!(A, z, F_1, F_2)
129+
A[1,2] = z*A[1,1]/(3im) + (F_1[1,1]-2F_2[1,1])/(3im)
130+
for j = 1:size(A,2)-2
131+
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)
132+
end
133+
A
134+
end
135+
136+
137+
"""
138+
logkernelsquare(z, n)
139+
140+
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``.
141+
The bottom right of the returned matrix is zero. For a square truncation compute `logkernelsquare(z,2n-1)[1:n,1:n]`.
142+
"""
143+
144+
function logkernelsquare(z, n)
145+
T = complex(float(eltype(z)))
146+
logkernelsquare!(zeros(T,n,n), z, zeros(T,n,n), zeros(T,n,n))
147+
end
148+
149+
150+
function logkernelsquare!(A::AbstractMatrix{T}, z, F_1, F_2) where T
151+
m,n = size(A)
152+
@assert m == n
153+
rec_rhs_1!(F_1, z)
154+
rec_rhs_2!(F_2, z)
155+
A[1,1] = L₀₀(z)
156+
logkernelsquare_populatefirstcolumn!(A, z, F_1, F_2)
157+
logkernelsquare_populatefirstrow!(A, z, F_1, F_2)
158+
159+
# F = F_1 # reuse the memory
160+
F = F_2 .- F_1
161+
162+
# 2nd row/column
163+
164+
for k = 1:m-2
165+
A[k+1,2] = im*(F[k+1,1] + (A[k,1] - A[k+2,1])/(2k+1))
166+
end
167+
168+
for j = 2:n-2
169+
A[2,j+1] = F[1,j+1] + im*(A[1,j+2] - A[1,j])/(2j+1)
170+
end
171+
172+
for ℓ = 1:((n-1)÷2-1)
173+
for k = ℓ+1:n-(ℓ+2)
174+
A[k+1,ℓ+2] = (2ℓ+1)*im*(F[k+1,ℓ+1] + (A[k,ℓ+1] - A[k+2,ℓ+1])/(2k+1)) + A[k+1,ℓ]
175+
end
176+
for j = ℓ+2:n-(ℓ+2)
177+
A[ℓ+2,j+1] = (2ℓ+1) * (F[ℓ+1,j+1] + im*(A[ℓ+1,j+2] - A[ℓ+1,j])/(2j+1)) + A[ℓ,j+1]
178+
end
179+
end
180+
A
181+
end

0 commit comments

Comments
 (0)