In [2]:
using ApproxFun, SingularIntegralEquations, RiemannHilbert, Plots, QuadGK, DualNumbers, ComplexPhasePortrait, LinearAlgebra
import ApproxFunBase: SequenceSpace, BasisFunctional, ℓ⁰, SpaceOperator, piece, pieces, npieces
import ApproxFunOrthogonalPolynomials: Recurrence
import RiemannHilbert: RiemannDual, logpart, finitepart, istieltjes, LogNumber
import SingularIntegralEquations: ⁺, ⁻
In [3]:
Γ = 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[3]:
We now define the jump:
In [4]:
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 [8]:
Φ = transpose(rhsolve(transpose(G), 800));
In [9]:
Φ(100.0)
Out[9]:
In [10]:
xx = range(-3.5,3.5,length=200)
Z = Φ[1,2].(xx' .+ im.*xx)
portrait(Z)
Out[10]:
In [11]:
ζ = exp(im*π/6)
Φ((ζ)⁺) - Φ((ζ)⁻)*G(ζ) |> norm
Out[11]:
To recover Painleve II, we actually want the limit:
In [12]:
10000.0Φ(10000.0)
Out[12]:
In [13]:
z = Fun(ℂ)
2(z*Φ[1,2])(Inf)
Out[13]:
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 [14]:
U = istieltjes(Φ)
plot(domain(U))
Out[14]:
In [15]:
I + stieltjes(U,10.0) == Φ(10.0)
Out[15]:
Thus we want to solve:
In [16]:
Φ((ζ)⁺) ≈ Φ((ζ)⁻)*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[16]:
We also have Plemelj that tells us ${\cal C}^+ -{\cal C^-} = I$, this changes (simplifies?) the equation to:
In [17]:
V(ζ) + cauchy(V,(ζ)⁻) - cauchy(V,(ζ)⁻)*G(ζ) ≈ G(ζ)-I
V(ζ) - cauchy(V,(ζ)⁻)*(G(ζ)-I) ≈ G(ζ)-I
Out[17]:
In [18]:
P = k -> Fun(Legendre(), [zeros(k);1])
xx = range(-2.,2.,length=200)
Z = xx' .+ im.*xx
portrait(stieltjes.(P(5), Z))
Out[18]:
In [19]:
portrait(Z.^(-6))
Out[19]:
Recall: the Jacobi operator is a tridiagonal operator that encodes multiplication by $x$:
In [20]:
x = Fun(Legendre())
J = Recurrence(Legendre())
f = Fun(exp, Legendre())
norm(J*f - x*f)
Out[20]:
The Jacobi operator also encodes the three-term recurrence:
In [21]:
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[21]:
The Stieltjes transform satisfies the exact same recurrence for $n \geq 1$:
In [22]:
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[22]:
But not for $n = 0$:
In [23]:
n = 0
z*Q(n,z) ≈ J[n+1,n+1]*Q(n,z) + J[n+2,n+1]*Q(n+1,z)
Out[23]:
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 [24]:
L = z -> [BasisFunctional(1);
BasisFunctional(2);
SpaceOperator(transpose(J - z*I)[2:end,:],ℓ⁰,ℓ⁰)]
L(0)
Out[24]:
In [25]:
μ₀ = stieltjes(P(0)*w,z)
μ₁ = stieltjes(P(1)*w,z)
n = 5
u = L(z)[1:n,1:n] \ [μ₀; μ₁; zeros(n-2)]
Out[25]:
This works great on or near the interval:
In [26]:
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[26]:
In [27]:
stieltjes(P(n-1),z)-u[end]
Out[27]:
Note this is much better than quadrature:
In [28]:
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 [29]:
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[29]:
In [30]:
stieltjes(P(n-1), z)
Out[30]:
But note that the true stieltjes is in $\ell^2$:
In [31]:
[stieltjes(P(k), z) for k=0:100] |> norm
Out[31]:
Thus we can solve it using one boundary condition:
In [32]:
T = z -> [BasisFunctional(1);
SpaceOperator(transpose(J - z*I)[2:end,:],ℓ⁰,ℓ⁰)]
T(z)
Out[32]:
In [38]:
z = 100.0
μ₀ = stieltjes(P(0)*w,z)
n = 100
u = T(z)[1:n,1:n] \ [μ₀;zeros(n-1)]
Out[38]:
In [39]:
norm(u - [stieltjes(P(k), z) for k=0:n-1])
Out[39]:
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 [40]:
f = Fun(exp, -1..1)
stieltjes(f, 1.000001)
Out[40]:
In [41]:
stieltjes(f, 1+0.000001im)
Out[41]:
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 [42]:
x = dual(1,2)
cos(x)
Out[42]:
In [43]:
cos(1), -2sin(1)
Out[43]:
In [44]:
cos(cos(x))
Out[44]:
In [45]:
cos(cos(1)), 2sin(cos(1))*sin(1)
Out[45]:
This is not numerical differentation: it is numerically stable.
In [46]:
# cos^n(x)
n = 100
ret = sin(x)
for k=1:n-1
ret = sin(ret)
end
ret
Out[46]:
In [47]:
dualpart(ret) - 2*0.0036986197855601063
Out[47]:
The direction can be complex as well:
In [48]:
cos(dual(1.0,im))
Out[48]:
In [49]:
sqrt(dual(-1.0+0im,im))
Out[49]:
In [50]:
sqrt(dual(-1.0-0im,-im))
Out[50]:
In [51]:
log(dual(0.0,1.0))
Out[51]:
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 [52]:
log(RiemannDual(0.0,1.0))
Out[52]:
In [53]:
l = log(RiemannDual(0.0,im))
Out[53]:
In [54]:
h = 0.00001
logpart(l)*log(h) + finitepart(l)
Out[54]:
In [55]:
l = stieltjes(f, RiemannDual(1.0,im))
l(h)
Out[55]:
In [56]:
stieltjes(f, 1.0 + im*h)
Out[56]:
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 [57]:
Γ = Segment(0,1) ∪ Segment(0, im)
plot(Γ;xlims=(-2,2),ylims=(-2,2))
Out[57]:
And here's a function that satisfies the zero-sum condition:
In [58]:
f = Fun(z -> angle(z) == 0 ? exp(z) : -cos(abs(z)), Γ)
f₁, f₂ = pieces(f)
first(f₁) + first(f₂)
Out[58]:
It's Stieltjes transform does not blow up near zero:
In [59]:
xx = range(-2,2,length=200)
Z = xx' .+ im*xx
portrait(stieltjes.(f,Z))
Out[59]:
Each component has a logarithmic singularity:
In [60]:
stieltjes(f₁, -0.000001)
Out[60]:
In [61]:
stieltjes(f₂, -0.000001)
Out[61]:
But these are cancelled when summed together:
In [62]:
stieltjes(f, -0.000001)
Out[62]:
This is easy to see with log numbers:
In [63]:
stieltjes(f₁, RiemannDual(0,-1))
Out[63]:
In [64]:
stieltjes(f₂, RiemannDual(0,-1))
Out[64]:
In [65]:
stieltjes(f, RiemannDual(0,-1))
Out[65]:
The value still depends on the angle:
In [66]:
finitepart(stieltjes(f,RiemannDual(0,1+im))),
finitepart(stieltjes(f,RiemannDual(0,-1+im)))
Out[66]:
In [67]:
G₁, G₂, G₃, G₄ = first.(pieces(G))
G₄*G₃*G₂*G₁
Out[67]:
The cyclic junction condition is a property satisfied by the right-hand side:
In [68]:
F = G-I
F₁, F₂, F₃, F₄ = first.(pieces(F))
F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄ |> norm
Out[68]:
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 [69]:
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[69]:
Thus combined with the cyclical jump condition, we satisfy the jumps:
In [70]:
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[70]:
In [71]:
norm( F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄)
Out[71]:
In [72]:
H = Hilbert() : ChebyshevWeight(Segment(-1,1),0)
Out[72]:
Thus $ I - (G-I)C^-$ can be represented on a single interval when $G$ has appropriate decay to 1:
In [73]:
G̃ = Fun(x->exp(-40x^2),ChebyshevWeight(Segment(-1,1),0))
G = Fun(x->1+exp(-40x^2),Chebyshev())
(I+G)/2 + im/2*G̃*H
Out[73]:
How to do this with junctions?
In [74]:
Fun(z -> stieltjes(P(15),z), Legendre(Segment(1+im,2+2im))) |>norm
Out[74]:
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.