In [3]:
using ApproxFun, SingularIntegralEquations, RiemannHilbert, Plots, QuadGK, DualNumbers, ComplexPhasePortrait
import ApproxFunBase: SequenceSpace, BasisFunctional, ℓ⁰, SpaceOperator, piece, pieces, npieces
import ApproxFunOrthogonalPolynomials: Recurrence
import RiemannHilbert: RiemannDual, logpart, finitepart, istieltjes, LogNumber
import SingularIntegralEquations: ⁺, ⁻
In [67]:
Γ = Segment(0, 2.5exp(im*π/6)) ∪
Segment(0, 2.5exp(5im*π/6)) ∪
Segment(0, 2.5exp(-5im*π/6)) ∪
Segment(0, 2.5exp(-im*π/6))
plot(Γ)
Out[67]:
We now define the jump:
In [71]:
s₁ = im
s₃ = -im
G = Fun( z -> if angle(z) ≈ π/6
[1 0;
s₁*exp(8im/3*z^3) 1]
elseif angle(z) ≈ 5π/6
[1 0;
s₃*exp(8im/3*z^3) 1]
elseif angle(z) ≈ -π/6
[1 -s₃*exp(-8im/3*z^3);
0 1]
elseif angle(z) ≈ -5π/6
[1 -s₁*exp(-8im/3*z^3);
0 1]
end
, Γ);
We can now solve the Riemann–Hilbert problem. The default is to solve $\Phi_+ = G \Phi_-$, so we transpose twice:
In [72]:
Φ = rhsolve(G.', 800).';
In [73]:
Φ(100.0)
Out[73]:
In [5]:
xx = linspace(-3.5,3.5,200)
Z = Φ[1,2].(xx' .+ im.*xx)
portrait(Z)
Out[5]:
In [78]:
ζ = exp(im*π/6)
Φ(ζ*⁺) - Φ(ζ*⁻)*G(ζ)|>norm
Out[78]:
To recover Painleve II, we actually want the limit:
In [81]:
10000.0Φ(10000.0)
Out[81]:
In [82]:
z = Fun(ℂ)
2(z*Φ[1,2])(Inf)
Out[82]:
The Stieltjes transform is defined via $$ {\cal S}_\Gamma f(z) = \int_\Gamma {f(\zeta) \over z - \zeta} {\rm d} \zeta $$ And the Cauchy transform is $$ {\cal C}_\Gamma f(z) = -{1 \over 2 \pi {\rm i}} {\cal S}_\Gamma(z) $$
The solution is given by $$ \Phi(z) = I + {\cal S}_\Gamma U(z) $$ for some $U$ defined on $\Gamma$:
In [85]:
U = istieltjes(Φ)
plot(domain(U))
Out[85]:
In [87]:
I + stieltjes(U,10.0) == Φ(10.0)
Out[87]:
Thus we want to solve:
In [88]:
Φ(ζ*⁺) ≈ Φ(ζ*⁻)*G(ζ)
I + stieltjes(U,ζ*⁺) ≈ (I + stieltjes(U,ζ*⁻))*G(ζ)
stieltjes(U,ζ*⁺) - stieltjes(U,ζ*⁻)*G(ζ) ≈ G(ζ)-I
V = U*(-2π*im)
cauchy(V,ζ*⁺) - cauchy(V,ζ*⁻)*G(ζ) ≈ G(ζ)-I
Out[88]:
We also have Plemelj that tells us ${\cal C}^+ -{\cal C^-} = I$, this changes (simplifies?) the equation to:
In [89]:
V(ζ) + cauchy(V,ζ*⁻) - cauchy(V,ζ*⁻)*G(ζ) ≈ G(ζ)-I
V(ζ) - cauchy(V,ζ*⁻)*(G(ζ)-I) ≈ G(ζ)-I
Out[89]:
In [90]:
P = k -> Fun(Legendre(), [zeros(k);1])
xx = linspace(-2.,2.,200)
Z = xx' .+ im.*xx
portrait(stieltjes.(P(5), Z))
Out[90]:
In [92]:
portrait(Z.^(-6))
Out[92]:
Recall: the Jacobi operator is a tridiagonal operator that encodes multiplication by $x$:
In [93]:
x = Fun(Legendre())
J = Recurrence(Legendre())
f = Fun(exp, Legendre())
norm(J*f - x*f)
Out[93]:
The Jacobi operator also encodes the three-term recurrence:
In [15]:
n=0
x*P(n) ≈ J[n+1,n+1]*P(n) + J[n+2,n+1]*P(n+1)
n=1
x*P(n) ≈ J[n,n+1]*P(n-1) + J[n+1,n+1]*P(n) + J[n+2,n+1]*P(n+1)
Out[15]:
The Stieltjes transform satisfies the exact same recurrence for $n \geq 1$:
In [94]:
z = 1.0+2.0im
w = 1 # in general, the weight corresponding to OP
Q = (n,z) -> sum(P(n)*w/(z-x)) # sum means integral
n = 1
z*Q(n,z) ≈ J[n,n+1]*Q(n-1,z) + J[n+1,n+1]*Q(n,z) + J[n+2,n+1]*Q(n+1,z)
Out[94]:
But not for $n = 0$:
In [17]:
n = 0
z*Q(n,z) ≈ J[n+1,n+1]*Q(n,z) + J[n+2,n+1]*Q(n+1,z)
Out[17]:
But if we know Q(0,z) and Q(1,z) in closed form (We do!), we can calculate Q via recurrence.
I like to represent this as a lower triangular operator:
In [18]:
L = z -> [BasisFunctional(1);
BasisFunctional(2);
SpaceOperator((J - z*I).'[2:end,:],ℓ⁰,ℓ⁰)]
L(0)
Out[18]:
In [19]:
μ₀ = stieltjes(P(0)*w,z)
μ₁ = stieltjes(P(1)*w,z)
n = 5
u = L(z)[1:n,1:n] \ [μ₀; μ₁; zeros(n-2)]
Out[19]:
This works great on or near the interval:
In [97]:
z = 0.1+eps()im
μ₀ = stieltjes(P(0)*w,z)
μ₁ = stieltjes(P(1)*w,z)
n = 100
u = L(z)[1:n,1:n] \ [μ₀; μ₁; zeros(n-2)]
Out[97]:
In [98]:
stieltjes(P(n-1),z)-u[end]
Out[98]:
Note this is much better than quadrature:
In [22]:
z=0.1+0.001im
@time μ=stieltjes(P(n-1),z)
@time sum(P(n-1)/(z-x))
@show abs(sum(P(n-1)/(z-x))-μ)
@time quadgk(x->P(n-1)(x)/(z-x),-1,1)
@show abs(quadgk(x->P(n-1)(x)/(z-x),-1,1)[1] - μ);
Unfortunately, this breaks down for large $z$:
In [23]:
z = 100.0
μ₀ = stieltjes(P(0)*w,z)
μ₁ = stieltjes(P(1)*w,z)
n = 100
u = L(z)[1:n,1:n] \ [μ₀; μ₁; zeros(n-2)]
Out[23]:
In [24]:
stieltjes(P(n-1), z)
Out[24]:
But note that the true stieltjes is in $\ell^2$:
In [25]:
[stieltjes(P(k), z) for k=0:100] |>norm
Out[25]:
Thus we can solve it using one boundary condition:
In [26]:
T = z -> [BasisFunctional(1);
SpaceOperator((J - z*I).'[2:end,:],ℓ⁰,ℓ⁰)]
T(z)
Out[26]:
In [105]:
z = 100.0
μ₀ = stieltjes(P(0)*w,z)
n = 100
u = T(z)[1:n,1:n] \ [μ₀;zeros(n-1)]
Out[105]:
In [106]:
norm(u - [stieltjes(P(k), z) for k=0:n-1])
Out[106]:
In practice, this can be determined adaptively using adaptive QR or F.W.J. Olver's algorithm (adaptive Gaussian elimination), but that's getting away.
Thus we have this algorithm:
In [108]:
f = Fun(exp, -1..1)
stieltjes(f, 1.000001)
Out[108]:
In [110]:
stieltjes(f, 1+0.000001im)
Out[110]:
The value has a "log part" which blows up, and a "finite part" that depends on the angle of approach. How do we represent this on a computer?
For inspiration, we use dual numbers:
In [31]:
x = dual(1,2)
cos(x)
Out[31]:
In [32]:
cos(1), -2sin(1)
Out[32]:
In [33]:
cos(cos(x))
Out[33]:
In [34]:
cos(cos(1)), 2sin(cos(1))*sin(1)
Out[34]:
This is not numerical differentation: it is numerically stable.
In [35]:
# cos^n(x)
n = 100
ret = sin(x)
for k=1:n-1
ret = sin(ret)
end
ret
Out[35]:
In [36]:
dualpart(ret) - 2*0.0036986197855601063
Out[36]:
The direction can be complex as well:
In [37]:
cos(dual(1.0,im))
Out[37]:
In [38]:
sqrt(dual(-1.0+0im,im))
Out[38]:
In [39]:
sqrt(dual(-1.0-0im,-im))
Out[39]:
In [40]:
log(dual(0.0,1.0))
Out[40]:
We use a dual-like number to annotate a number with a direction. When evaluated on a logarithmic singularity, it returns a "Log number", which represents local logarithmic behaviour of a function $f(z)$ by two numbers $a$ and $b$ so that $$ f(c + hd) \sim a \log h + b $$
In [111]:
log(RiemannDual(0.0,1.0))
Out[111]:
In [42]:
l = log(RiemannDual(0.0,im))
Out[42]:
In [43]:
h = 0.00001
logpart(l)*log(h) + finitepart(l)
Out[43]:
In [44]:
l = stieltjes(f, RiemannDual(1.0,im))
l(h)
Out[44]:
In [45]:
stieltjes(f, 1.0 + im*h)
Out[45]:
When we have contours meeting at a junction, it's possible for the logarithmic singularities to cancel. This is the case if they satisfy the "zero-sum condition": that the values of the limits sum to zero.
Here's a simple example of two contours who have 0 as a joint junction point:
In [46]:
Γ = Segment(0,1) ∪ Segment(0, im)
plot(Γ;xlims=(-2,2),ylims=(-2,2))
Out[46]:
And here's a function that satisfies the zero-sum condition:
In [47]:
f = Fun(z -> angle(z) == 0 ? exp(z) : -cos(abs(z)), Γ)
f₁, f₂ = pieces(f)
first(f₁) + first(f₂)
Out[47]:
It's Stieltjes transform does not blow up near zero:
In [48]:
xx = linspace(-2,2,200)
Z = xx' .+ im*xx
portrait(stieltjes.(f,Z))
Out[48]:
Each component has a logarithmic singularity:
In [49]:
stieltjes(f₁, -0.000001)
Out[49]:
In [50]:
stieltjes(f₂, -0.000001)
Out[50]:
But these are cancelled when summed together:
In [51]:
stieltjes(f, -0.000001)
Out[51]:
This is easy to see with log numbers:
In [52]:
stieltjes(f₁, RiemannDual(0,-1))
Out[52]:
In [53]:
stieltjes(f₂, RiemannDual(0,-1))
Out[53]:
In [54]:
stieltjes(f, RiemannDual(0,-1))
Out[54]:
The value still depends on the angle:
In [55]:
finitepart(stieltjes(f,RiemannDual(0,1+im))),
finitepart(stieltjes(f,RiemannDual(0,-1+im)))
Out[55]:
In [56]:
G₁, G₂, G₃, G₄ = first.(pieces(G))
G₄*G₃*G₂*G₁
Out[56]:
The cyclic junction condition is a property satisfied by the right-hand side:
In [57]:
F = G-I
F₁, F₂, F₃, F₄ = first.(pieces(F))
F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄ |> norm
Out[57]:
This is due to a telescoping sum argument. In general, any function U that satisfies the zero sum condition has the property that $(I - (G-I)C^-) U $ satisfies the cyclic junction condition.
Zero-sum condition also gets us continuity in each sector:
In [58]:
z₁, z₂, z₃, z₄ = exp(π*im/6), exp(5π*im/6), exp(-5π*im/6), exp(-π*im/6)
U = istieltjes(Φ)
U₁, U₂, U₃, U₄ = pieces(U)
C₁₊, C₂₊, C₃₊, C₄₊ = finitepart.(
stieltjes(U₁,RiemannDual(0,z₁)⁺) +
stieltjes(U₂,RiemannDual(0,z₁)) +
stieltjes(U₃,RiemannDual(0,z₁)) +
stieltjes(U₄,RiemannDual(0,z₁))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₂)) +
stieltjes(U₂,RiemannDual(0,z₂)⁺) +
stieltjes(U₃,RiemannDual(0,z₂)) +
stieltjes(U₄,RiemannDual(0,z₂))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₃)) +
stieltjes(U₂,RiemannDual(0,z₃)) +
stieltjes(U₃,RiemannDual(0,z₃)⁺) +
stieltjes(U₄,RiemannDual(0,z₃))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₄)) +
stieltjes(U₂,RiemannDual(0,z₄)) +
stieltjes(U₃,RiemannDual(0,z₄)) +
stieltjes(U₄,RiemannDual(0,z₄)⁺)
)
C₁₋, C₂₋, C₃₋, C₄₋ = finitepart.(
stieltjes(U₁,RiemannDual(0,z₁)⁻) +
stieltjes(U₂,RiemannDual(0,z₁)) +
stieltjes(U₃,RiemannDual(0,z₁)) +
stieltjes(U₄,RiemannDual(0,z₁))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₂)) +
stieltjes(U₂,RiemannDual(0,z₂)⁻) +
stieltjes(U₃,RiemannDual(0,z₂)) +
stieltjes(U₄,RiemannDual(0,z₂))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₃)) +
stieltjes(U₂,RiemannDual(0,z₃)) +
stieltjes(U₃,RiemannDual(0,z₃)⁻) +
stieltjes(U₄,RiemannDual(0,z₃))
) ,
finitepart.(
stieltjes(U₁,RiemannDual(0,z₄)) +
stieltjes(U₂,RiemannDual(0,z₄)) +
stieltjes(U₃,RiemannDual(0,z₄)) +
stieltjes(U₄,RiemannDual(0,z₄)⁻)
) ;
norm(C₁₊ - C₂₋)
Out[58]:
Thus combined with the cyclical jump condition, we satisfy the jumps:
In [59]:
C₄₊ ≈ C₄₋*G₄ + F₄
C₄₊ ≈ C₃₊*G₄ + F₄
C₄₊ ≈ C₃₋*G₃*G₄ + F₃*G₄ + F₄
C₄₊ ≈ C₂₊*G₃*G₄ + F₃*G₄ + F₄
C₄₊ ≈ C₁₊*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄
C₄₊ ≈ C₁₋*G₁*G₂*G₃*G₄ + F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄
Out[59]:
In [60]:
norm( F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄)
Out[60]:
In [61]:
H = Hilbert() : ChebyshevWeight()
Out[61]:
Thus $ I - (G-I)C^-$ can be represented on a single interval when $G$ has appropriate decay to 1:
In [62]:
G̃ = Fun(x->exp(-40x^2),ChebyshevWeight())
G = Fun(x->1+exp(-40x^2),Chebyshev())
(I+G)/2 + im/2*G̃*H
Out[62]:
How to do this with junctions?
In [63]:
Fun(z -> stieltjes(P(15),z), Legendre(Segment(1+im,2+2im))) |>norm
Out[63]:
In [64]:
K = LowRankFun((z,ζ) -> 1/(z-ζ), Legendre()*Legendre(Segment(1+im,2+2im)))
rank(K)
Out[64]:
This is exactly the sort of structure that makes Fast Multipole Method / Hierarchical matrices fast. We can exploit this by plugging in RiemannHilbert.jl into HierarchicalMatrices.jl.