# Riemann–Hilbert problems



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



# 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 [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]:

-2

-1

0

1

2

-1.0

-0.5

0.0

0.5

1.0

y1



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]:

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 [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]:

8.263986678835547e-16



To recover Painleve II, we actually want the limit:



In [12]:

10000.0Φ(10000.0)




Out[12]:

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 [13]:

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




Out[13]:

-0.36706155154807807 - 6.938893903907228e-17im



## 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 [14]:

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




Out[14]:

-2

-1

0

1

2

-1.0

-0.5

0.0

0.5

1.0

y1




In [15]:

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




Out[15]:

true



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]:

true



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]:

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 [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]:

0.0



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]:

true



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]:

true



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]:

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 [24]:

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

L(0)




Out[24]:

InterlaceOperator : SequenceSpace() → 3-element ArraySpace:
Space[ConstantSpace, ConstantSpace, SequenceSpace()]
1.0                 0.0                  …  0.0                 ⋯
0.0                 1.0                     0.0                 ⋱
0.3333333333333333  0.0                     0.0                 ⋱
0.0                 0.39999999999999997     0.0                 ⋱
0.0                 0.0                     0.0                 ⋱
0.0                 0.0                  …  0.0                 ⋱
0.0                 0.0                     0.0                 ⋱
0.0                 0.0                     0.0                 ⋱
0.0                 0.0                     0.0                 ⋱
0.0                 0.0                     0.5294117647058824  ⋱
⋮                   ⋱                   …   ⋱                  ⋱




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]:

5-element Array{Complex{Float64},1}:
0.34657359027997253 - 0.7853981633974483im
-0.08263008292513074 - 0.09225098283750298im
-0.02047897101517343 + 0.006432358667077497im
-0.0004860919654602121 + 0.0039580162862197284im
0.000655510320055656 + 0.0004009376214656593im



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]:

100-element Array{Complex{Float64},1}:
0.20067069546215124 - 3.1415926535897927im
-1.9799329304537843 - 0.3141592653589794im
-0.3973252872991432 + 1.5236724369910486im
1.2537344057526651 + 0.4633849164044942im
0.5173974864810735 - 1.0616619673724992im
-0.9098559770355384 - 0.561807087250645im
-0.597971501190743 + 0.7817203401477975im
0.6688247015236087 + 0.6267255665280007im
0.6486296950775765 - 0.5664942539053223im
-0.4719919033952208 - 0.66409386265145im
-0.6734451872149104 + 0.38366699461101433im
0.3005167400728086 + 0.6769672104725115im
0.6749237967942893 - 0.22194269638619826im
⋮
0.14184941367271484 + 0.22634421744973532im
0.23805329322050453 - 0.11784046844462849im
-0.09292715409138425 - 0.24726644153539204im
-0.25392063603186216 + 0.06736395129918642im
0.04140894981448789 + 0.25797833139669324im
0.25942757082257345 - 0.015321337512756154im
0.01064109831526331 - 0.25828184927001063im
-0.25457973578066756 - 0.0362244339474118im
-0.06118101347250797 + 0.2483842603359941im
0.23978207293736778 + 0.0852717719941878im
0.1082684564195096 - 0.22888238468066074im
-0.2158157025898757 - 0.12995572357211532im




In [27]:

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




Out[27]:

2.7755575615628914e-17 + 1.1102230246251565e-16im



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




0.000079 seconds (31 allocations: 6.375 KiB)
0.067869 seconds (819 allocations: 24.827 MiB, 26.75% gc time)
abs(sum(P(n - 1) / (z - x)) - μ) = 4.392606376355813e-14
1.357889 seconds (2.46 M allocations: 128.622 MiB, 5.79% gc time)
#= In[28]:6 =#
(P(n - 1))(x) / (z - x)
end), -1, 1))[1] - μ) = 4.3885418357208767e-16



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]:

100-element Array{Float64,1}:
0.020000666706669522
6.667066695240317e-5
2.6668952571412785e-7
1.142984085868605e-9
5.0707414099784735e-12
-1.6538148987588797e-12
-3.074250159474433e-10
-5.709179969175481e-8
-1.0704443445315075e-5
-0.002021899680293121
-0.38415130525659225
-73.33613836745825
-14055.741048399677
⋮
-3.9496900117655246e178
-7.854805097786174e180
-1.5621944002918002e183
-3.1071441442662403e185
-6.180360461213906e187
-1.2293958060198338e190
-2.445651787132766e192
-4.865438225763328e194
-9.679952785716854e196
-1.925963071583903e199
-3.8321776469863156e201
-7.625455777961642e203




In [30]:

stieltjes(P(n-1), z)




Out[30]:

0.0



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



In [31]:

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




Out[31]:

0.020000777828880434



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]:

InterlaceOperator : SequenceSpace() → 2-element ArraySpace:
Space[ConstantSpace, SequenceSpace()]
1.0                    0.0                  …     0.0                 ⋯
0.3333333333333333  -100.0                        0.0                 ⋱
0.0                    0.39999999999999997        0.0                 ⋱
0.0                    0.0                        0.0                 ⋱
0.0                    0.0                        0.0                 ⋱
0.0                    0.0                  …     0.0                 ⋱
0.0                    0.0                        0.0                 ⋱
0.0                    0.0                        0.0                 ⋱
0.0                    0.0                        0.5294117647058824  ⋱
0.0                    0.0                     -100.0                 ⋱
⋮                      ⋱                   …      ⋱                  ⋱




In [38]:

z = 100.0

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

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




Out[38]:

100-element Array{Float64,1}:
0.020000666706669522
6.667066695240317e-5
2.666895257144473e-7
1.1429841391064577e-9
5.080057794657254e-12
2.3091753139353408e-14
1.065800004129641e-16
4.973858189823564e-19
2.3406978482382525e-21
1.1087794019770231e-23
5.280034191637557e-26
2.525296987630347e-28
1.2121729040401978e-30
⋮
6.09287145569266e-206
3.0294922753741863e-208
1.5064150369499735e-210
7.491103591762882e-213
3.7253987099462816e-215
1.852784717942273e-217
9.215138511587441e-220
4.583560451533549e-222
2.2799627188325937e-224
1.1341636567013033e-226
5.6421734618294046e-229
2.806910415683975e-231




In [39]:

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




Out[39]:

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 [40]:

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




Out[40]:

35.85252983289078




In [41]:

stieltjes(f, 1+0.000001im)




Out[41]:

35.852495347781634 - 4.269828356622018im



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]:

0.5403023058681398 - 1.682941969615793ɛ




In [43]:

cos(1), -2sin(1)




Out[43]:

(0.5403023058681398, -1.682941969615793)




In [44]:

cos(cos(x))




Out[44]:

0.8575532158463934 + 0.8656973695406469ɛ




In [45]:

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




Out[45]:

(0.8575532158463934, 0.8656973695406469)



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]:

0.16885248872798114 + 0.007397239571120209ɛ




In [47]:

dualpart(ret) - 2*0.0036986197855601063




Out[47]:

-3.469446951953614e-18



The direction can be complex as well:



In [48]:

cos(dual(1.0,im))




Out[48]:

0.5403023058681398 - 0.0im + 0.0ɛ - 0.8414709848078965imɛ




In [49]:

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




Out[49]:

0.0 + 1.0im + 0.5ɛ + 0.0imɛ




In [50]:

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




Out[50]:

0.0 - 1.0im + 0.5ɛ + 0.0imɛ




In [51]:

log(dual(0.0,1.0))




Out[51]:

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 [52]:

log(RiemannDual(0.0,1.0))




Out[52]:

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




In [53]:

l = log(RiemannDual(0.0,im))




Out[53]:

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




In [54]:

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




Out[54]:

-11.512925464970229 + 1.5707963267948966im




In [55]:

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

l(h)




Out[55]:

29.593415861168136 - 4.269867111336779im




In [56]:

stieltjes(f, 1.0 + im*h)




Out[56]:

29.59345855815111 - 4.269542154749188im



## 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 [57]:

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

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




Out[57]:

-2

-1

0

1

2

-2

-1

0

1

2

y1



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]:

-5.551115123125783e-16 + 0.0im



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]:

-15.133399294132499 + 0.0im




In [61]:

stieltjes(f₂, -0.000001)




Out[61]:

13.575698815971109 + 1.5707948403214074im



But these are cancelled when summed together:



In [62]:

stieltjes(f, -0.000001)




Out[62]:

-1.5577004781613901 + 1.5707948403214074im



This is easy to see with log numbers:



In [63]:

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




Out[63]:

(1.0000000000000007 + 0.0im)log ε + -1.3179021514543987 + 0.0im




In [64]:

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




Out[64]:

(-1.0000000000000013 + 0.0im)log ε + -0.23981174200057068 + 1.5707963267948988im




In [65]:

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




Out[65]:

(-6.661338147750939e-16 + 0.0im)log ε + -1.5577138934549695 + 1.5707963267948988im



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]:

(-1.5577138934549697 - 4.712388980384695im, -1.5577138934549697 + 1.5707963267948983im)



## 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 [67]:

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

G₄*G₃*G₂*G₁




Out[67]:

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 [68]:

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

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




Out[68]:

9.675222525179979e-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 [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]:

6.214465357936134e-16



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]:

true




In [71]:

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




Out[71]:

9.675222525179979e-18



## 2017 v 2023 (Future directions.)

President Ivanka Trump.... 😰

#### Banded operators

Chebyshev singularities lead to banded operators:



In [72]:

H = Hilbert() : ChebyshevWeight(Segment(-1,1),0)




Out[72]:

HilbertWrapper : (1-x^2)^-0.5[Chebyshev(the segment [-1.0,1.0])] → Ultraspherical(1,the segment [-1.0,1.0])
0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  ⋱
0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  ⋱
0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  ⋱
0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  ⋱
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  ⋱
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  ⋱
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.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 [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]:

PlusOperator : (1-x^2)^-0.5[Chebyshev(the segment [-1.0,1.0])] → (1-x^2)^-0.5[Ultraspherical(1,the segment [-1.0,1.0])]
1.0886432670340571 + 0.0im  …  ⋯
0.0 + 0.0im     ⋱
-0.08426795592489292 + 0.0im     ⋱
0.0 + 0.0im     ⋱
0.0761649869582432 + 0.0im     ⋱
0.0 + 0.0im  …  ⋱
-0.06546994432074481 + 0.0im     ⋱
0.0 + 0.0im     ⋱
0.053542023442444714 + 0.0im     ⋱
0.0 + 0.0im     ⋱
⋮       …  ⋱



How to do this with junctions?

#### Hierarchical structure

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



In [74]:

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




Out[74]:

8.244707787129414e-9




In [64]:

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

rank(K)




Out[64]:

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.