# Riemann–Hilbert problems-checkpoint



In :

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: ⁺, ⁻



# Riemann–Hilbert problems

Sheehan Olver
Imperial College

## RiemannHilbert.jl

An in-development package for solving Riemann–Hilbert problems, that builds on ApproxFun.jl and SingularIntegralEquations.jl.

Here's an example of calculating the Hastings–McLeod solution. This has the following contour:



In :

Γ = 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:

-2

-1

0

1

2

-1.0

-0.5

0.0

0.5

1.0

y1



We now define the jump:



In :

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 :

Φ = rhsolve(G.', 800).';




In :

Φ(100.0)




Out:

2×2 Array{Complex{Float64},2}:
1.0-0.000345453im  -0.00183531+6.75051e-6im
-0.00183531-6.75051e-6im           1.0+0.000345453im




In :

xx = linspace(-3.5,3.5,200)
Z = Φ[1,2].(xx' .+ im.*xx)
portrait(Z)




Out:




In :

ζ = exp(im*π/6)
Φ(ζ*⁺) - Φ(ζ*⁻)*G(ζ)|>norm




Out:

4.4191721953081e-16



To recover Painleve II, we actually want the limit:



In :

10000.0Φ(10000.0)




Out:

2×2 Array{Complex{Float64},2}:
10000.0-0.0345457im   -0.183531+6.75028e-6im
-0.183531-6.75028e-6im    10000.0+0.0345457im




In :

z = Fun(ℂ)
2(z*Φ[1,2])(Inf)




Out:

-0.36706155154807807 + 2.220446049250313e-16im



## Stieltjes transforms, Cauchy transforms, Riemann–Hilbert problems and singular integral equations

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 :

U = istieltjes(Φ)
plot(domain(U))




Out:

-2

-1

0

1

2

-1.0

-0.5

0.0

0.5

1.0

y1




In :

I + stieltjes(U,10.0) == Φ(10.0)




Out:

true



Thus we want to solve:



In :

Φ(ζ*⁺) ≈ Φ(ζ*⁻)*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:

true



We also have Plemelj that tells us ${\cal C}^+ -{\cal C^-} = I$, this changes (simplifies?) the equation to:



In :

V(ζ) + cauchy(V,ζ*⁻) - cauchy(V,ζ*⁻)*G(ζ) ≈ G(ζ)-I
V(ζ)  - cauchy(V,ζ*⁻)*(G(ζ)-I) ≈ G(ζ)-I




Out:

true



## Stieltjes transforms and OPs

The Stieltjes transform of weighted orthogonal polynomials satisfy the same recurrence as OPs themselves.

We need to calculate Stieltjes transforms of OPs:



In :

P = k -> Fun(Legendre(), [zeros(k);1])

xx = linspace(-2.,2.,200)
Z = xx' .+ im.*xx

portrait(stieltjes.(P(5), Z))




Out:




In :

portrait(Z.^(-6))




Out:



Recall: the Jacobi operator is a tridiagonal operator that encodes multiplication by $x$:



In :

x = Fun(Legendre())
J = Recurrence(Legendre())

f = Fun(exp, Legendre())
norm(J*f - x*f)




Out:

0.0



The Jacobi operator also encodes the three-term recurrence:



In :

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:

true



The Stieltjes transform satisfies the exact same recurrence for $n \geq 1$:



In :

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:

true



But not for $n = 0$:



In :

n = 0
z*Q(n,z) ≈ J[n+1,n+1]*Q(n,z) + J[n+2,n+1]*Q(n+1,z)




Out:

false



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 :

L = z -> [BasisFunctional(1);
BasisFunctional(2);
SpaceOperator((J - z*I).'[2:end,:],ℓ⁰,ℓ⁰)]

L(0)




Out:

InterlaceOperator:ApproxFun.SequenceSpace()→3-element ArraySpace:
ConstantSpace
ConstantSpace
ApproxFun.SequenceSpace()
1.0       0.0  0.0       0.0       …
0.0       1.0  0.0       0.0
0.333333  0.0  0.666667  0.0
0.4  0.0       0.6          0.0
0.428571  0.0          0.0       0.0
0.444444  …  0.0       0.0       0.0
0.0       0.0       0.0       ⋱
0.538462  0.0       0.0       ⋱
0.0       0.533333  0.0       ⋱
0.470588  0.0       0.529412  ⋱
…             ⋱         ⋱        ⋱




In :

μ₀ = stieltjes(P(0)*w,z)
μ₁ = stieltjes(P(1)*w,z)

n = 5
u = L(z)[1:n,1:n] \ [μ₀; μ₁; zeros(n-2)]




Out:

5-element Array{Complex{Float64},1}:
0.346574-0.785398im
-0.0826301-0.092251im
-0.020479+0.00643236im
-0.000486092+0.00395802im
0.00065551+0.000400938im



This works great on or near the interval:



In :

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:

100-element Array{Complex{Float64},1}:
0.200671-3.14159im
-1.97993-0.314159im
-0.397325+1.52367im
1.25373+0.463385im
0.517397-1.06166im
-0.909856-0.561807im
-0.597972+0.78172im
0.668825+0.626726im
0.64863-0.566494im
-0.471992-0.664094im
-0.673445+0.383667im
0.300517+0.676967im
0.674924-0.221943im
⋮
0.141849+0.226344im
0.238053-0.11784im
-0.0929272-0.247266im
-0.253921+0.067364im
0.0414089+0.257978im
0.259428-0.0153213im
0.0106411-0.258282im
-0.25458-0.0362244im
-0.061181+0.248384im
0.239782+0.0852718im
0.108268-0.228882im
-0.215816-0.129956im




In :

stieltjes(P(n-1),z)-u[end]




Out:

2.7755575615628914e-17 + 1.1102230246251565e-16im



Note this is much better than quadrature:



In :

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))-μ)




0.000172 seconds (165 allocations: 14.670 KiB)
0.077479 seconds (457.23 k allocations: 56.141 MiB, 33.08% gc time)
abs(sum(P(n - 1) / (z - x)) - μ) = 2.204983093081963e-14
0.938067 seconds (1.00 M allocations: 54.327 MiB, 5.27% gc time)
(P(n - 1))(x) / (z - x)
end), -1, 1)) - μ) = 4.3885418357208767e-16



Unfortunately, this breaks down for large $z$:



In :

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:

100-element Array{Float64,1}:
0.0200007
6.66707e-5
2.6669e-7
1.14298e-9
5.0955e-12
2.80279e-12
5.09598e-10
9.46372e-8
1.7744e-5
0.00335157
0.636782
121.564
23299.2
⋮
6.54713e178
1.30204e181
2.58954e183
5.1505e185
1.02448e188
2.03789e190
4.05399e192
8.06511e194
1.60458e197
3.19254e199
6.35234e201
1.26402e204




In :

stieltjes(P(n-1), z)




Out:

0.0



But note that the true stieltjes is in $\ell^2$:



In :

[stieltjes(P(k), z) for k=0:100] |>norm




Out:

0.020000777828880434



Thus we can solve it using one boundary condition:



In :

T = z -> [BasisFunctional(1);
SpaceOperator((J - z*I).'[2:end,:],ℓ⁰,ℓ⁰)]
T(z)




Out:

InterlaceOperator:ApproxFun.SequenceSpace()→2-element ArraySpace:
ConstantSpace
ApproxFun.SequenceSpace()
1.0          0.0     0.0          0.0       …
0.333333  -100.0     0.666667     0.0
0.4  -100.0          0.6
0.428571  -100.0
0.444444
…     0.0
0.0          0.0
0.533333     0.0       ⋱
-100.0          0.529412  ⋱
0.473684  -100.0       ⋱
…                   ⋱        ⋱




In :

z = 100.0

μ₀ = stieltjes(P(0)*w,z)

n = 100
u = T(z)[1:n,1:n] \ [μ₀;zeros(n-1)]




Out:

100-element Array{Float64,1}:
0.0200007
6.66707e-5
2.6669e-7
1.14298e-9
5.08006e-12
2.30918e-14
1.0658e-16
4.97386e-19
2.3407e-21
1.10878e-23
5.28003e-26
2.5253e-28
1.21217e-30
⋮
6.09287e-206
3.02949e-208
1.50642e-210
7.4911e-213
3.7254e-215
1.85278e-217
9.21514e-220
4.58356e-222
2.27996e-224
1.13416e-226
5.64217e-229
2.80691e-231




In :

norm(u - [stieltjes(P(k), z) for k=0:n-1])




Out:

5.739713219974347e-23



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:

1. If $z$ is close to $[-1,1]$ (where closed depends on $n$), solve a lower triangular system using forward recurrence.
2. If $z$ is far from $[-1,1]$, solve a tridiagonal system using adaptive QR.

## Dual numbers and "Log numbers"

We want to capture the behaviour of the singularities to set up a collocation system:



In :

f = Fun(exp, -1..1)
stieltjes(f, 1.000001)




Out:

35.85252983289078




In :

stieltjes(f, 1+0.000001im)




Out:

35.85249534778163 - 4.269828356622017im



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 :

x = dual(1,2)
cos(x)




Out:

0.5403023058681398 - 1.682941969615793ɛ




In :

cos(1), -2sin(1)




Out:

(0.5403023058681398, -1.682941969615793)




In :

cos(cos(x))




Out:

0.8575532158463934 + 0.8656973695406469ɛ




In :

cos(cos(1)), 2sin(cos(1))*sin(1)




Out:

(0.8575532158463934, 0.8656973695406469)



This is not numerical differentation: it is numerically stable.



In :

# cos^n(x)
n = 100
ret = sin(x)
for k=1:n-1
ret = sin(ret)
end
ret




Out:

0.16885248872798114 + 0.007397239571120209ɛ




In :

dualpart(ret) - 2*0.0036986197855601063




Out:

-3.469446951953614e-18



The direction can be complex as well:



In :

cos(dual(1.0,im))




Out:

0.5403023058681398 - 0.0im + 0.0ɛ - 0.8414709848078965imɛ




In :

sqrt(dual(-1.0+0im,im))




Out:

0.0 + 1.0im + 0.5ɛ + 0.0imɛ




In :

sqrt(dual(-1.0-0im,-im))




Out:

0.0 - 1.0im + 0.5ɛ + 0.0imɛ




In :

log(dual(0.0,1.0))




Out:

Dual{Float64}(-Inf,Inf)



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 :

log(RiemannDual(0.0,1.0))




Out:

(1.0 + 0.0im)log ε + 0.0 + 0.0im




In :

l = log(RiemannDual(0.0,im))




Out:

(1.0 + 0.0im)log ε + 0.0 + 1.5707963267948966im




In :

h = 0.00001
logpart(l)*log(h) + finitepart(l)




Out:

-11.512925464970229 + 1.5707963267948966im




In :

l = stieltjes(f, RiemannDual(1.0,im))

l(h)




Out:

29.593415861168136 - 4.269867111336779im




In :

stieltjes(f, 1.0 + im*h)




Out:

29.593458558151106 - 4.26954215474919im



## Zero-sum condition

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 :

Γ = Segment(0,1) ∪ Segment(0, im)

plot(Γ;xlims=(-2,2),ylims=(-2,2))




Out:

-2

-1

0

1

2

-2

-1

0

1

2

y1



And here's a function that satisfies the zero-sum condition:



In :

f = Fun(z -> angle(z) == 0 ? exp(z) : -cos(abs(z)), Γ)
f₁, f₂ = pieces(f)
first(f₁) + first(f₂)




Out:

-4.440892098500626e-16 + 0.0im



It's Stieltjes transform does not blow up near zero:



In :

xx = linspace(-2,2,200)
Z = xx' .+ im*xx
portrait(stieltjes.(f,Z))




Out:



Each component has a logarithmic singularity:



In :

stieltjes(f₁, -0.000001)




Out:

-15.133399294132497 + 0.0im




In :

stieltjes(f₂, -0.000001)




Out:

13.575698815971107 + 1.5707948403214076im



But these are cancelled when summed together:



In :

stieltjes(f, -0.000001)




Out:

-1.5577004781613901 + 1.5707948403214076im



This is easy to see with log numbers:



In :

stieltjes(f₁, RiemannDual(0,-1))




Out:

(1.0000000000000009 + 0.0im)log ε + -1.3179021514543992 + 0.0im




In :

stieltjes(f₂, RiemannDual(0,-1))




Out:

(-1.000000000000001 + 0.0im)log ε + -0.23981174200057237 + 1.5707963267948983im




In :

stieltjes(f, RiemannDual(0,-1))




Out:

(-2.220446049250313e-16 + 0.0im)log ε + -1.5577138934549715 + 1.5707963267948983im



The value still depends on the angle:



In :

finitepart(stieltjes(f,RiemannDual(0,1+im))),
finitepart(stieltjes(f,RiemannDual(0,-1+im)))




Out:

(-1.5577138934549715 - 4.712388980384694im, -1.5577138934549715 + 1.5707963267948992im)



## Product condition and cyclic junction condition

Now consider a Riemann–Hilbert problem which satisfies the product condition. That is, the product of the jumps is I:



In :

G₁, G₂, G₃, G₄ = first.(pieces(G))

G₄*G₃*G₂*G₁




Out:

2×2 Array{Complex{Float64},2}:
1.0+0.0im  6.84142e-18+0.0im
-6.84142e-18+0.0im          1.0+0.0im



The cyclic junction condition is a property satisfied by the right-hand side:



In :

F = G-I
F₁, F₂, F₃, F₄ = first.(pieces(F))

F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄   |> norm




Out:

6.841415457043595e-18



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 :

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:

6.483657258559754e-16



Thus combined with the cyclical jump condition, we satisfy the jumps:



In :

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:

true




In :

norm( F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄)




Out:

6.841415457043595e-18



## 2017 v 2023 (Future directions.)

President Ivanka Trump.... 😰

#### Banded operators

Chebyshev singularities lead to banded operators:



In :

H = Hilbert() : ChebyshevWeight()




Out:

HilbertWrapper:(1-x^2)^-0.5[Chebyshev(【-1.0,1.0】)]→Ultraspherical(1,【-1.0,1.0】)
0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  1.0
0.0  0.0  ⋱
⋱   ⋱



Thus $I - (G-I)C^-$ can be represented on a single interval when $G$ has appropriate decay to 1:



In :

G̃ = Fun(x->exp(-40x^2),ChebyshevWeight())
G = Fun(x->1+exp(-40x^2),Chebyshev())
(I+G)/2 + im/2*G̃*H




Out:

PlusOperator:(1-x^2)^-0.5[Chebyshev(【-1.0,1.0】)]→(1-x^2)^-0.5[Ultraspherical(1,【-1.0,1.0】)]
1.08864+0.0im          0.0+0.0880911im  …         0.0+0.0535527im  ⋱
0.0+0.0im     0.502188+0.0im            0.0118945+0.0im        ⋱
-0.084268+0.0im          0.0-0.0837981im            0.0-0.0535835im  ⋱
0.0+0.0im  -0.00405148+0.0im           -0.0226296+0.0im        ⋱
0.076165+0.0im          0.0+0.0758392im            0.0+0.0533786im  ⋱
0.0+0.0im   0.00534752+0.0im        …   0.0312126+0.0im        ⋱
-0.0654699+0.0im          0.0-0.0653164im            0.0-0.052508im   ⋱
0.0+0.0im  -0.00596396+0.0im            -0.536959+0.0im        ⋱
0.053542+0.0im          0.0+0.0535527im            0.0+0.0505444im  ⋱
0.0+0.0im   0.00593056+0.0im             0.539582+0.0im        ⋱
⋱                  ⋱            …             ⋱            ⋱



How to do this with junctions?

#### Hierarchical structure

Stieltjes transforms have low rank structure when evaluated off the interval:



In :

Fun(z -> stieltjes(P(15),z), Legendre(Segment(1+im,2+2im))) |>norm




Out:

8.244707787129415e-9




In :

K = LowRankFun((z,ζ) -> 1/(z-ζ), Legendre()*Legendre(Segment(1+im,2+2im)))

rank(K)




Out:

12



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.