Riemann–Hilbert problems-checkpoint



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

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 [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]:
-2 -1 0 1 2 -1.0 -0.5 0.0 0.5 1.0 y1

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

To recover Painleve II, we actually want the limit:


In [81]:
10000.0Φ(10000.0)


Out[81]:
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 [82]:
z = Fun()
2(z*Φ[1,2])(Inf)


Out[82]:
-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 [85]:
U = istieltjes(Φ)
plot(domain(U))


Out[85]:
-2 -1 0 1 2 -1.0 -0.5 0.0 0.5 1.0 y1

In [87]:
I + stieltjes(U,10.0) == Φ(10.0)


Out[87]:
true

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

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

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

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

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]:
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 [18]:
L = z -> [BasisFunctional(1);
     BasisFunctional(2);
    SpaceOperator((J - z*I).'[2:end,:],ℓ⁰,ℓ⁰)]
        
L(0)


Out[18]:
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 [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]:
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 [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]:
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 [98]:
stieltjes(P(n-1),z)-u[end]


Out[98]:
2.7755575615628914e-17 + 1.1102230246251565e-16im

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


  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)
abs((quadgk((x->begin  # In[22], line 6:
                    (P(n - 1))(x) / (z - x)
                end), -1, 1))[1] - μ) = 4.3885418357208767e-16

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]:
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 [24]:
stieltjes(P(n-1), z)


Out[24]:
0.0

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


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


Out[25]:
0.020000777828880434

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]:
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 [105]:
z = 100.0

μ₀ = stieltjes(P(0)*w,z)        
        
n = 100
u = T(z)[1:n,1:n] \ [μ₀;zeros(n-1)]


Out[105]:
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 [106]:
norm(u - [stieltjes(P(k), z) for k=0:n-1])


Out[106]:
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 [108]:
f = Fun(exp, -1..1)
stieltjes(f, 1.000001)


Out[108]:
35.85252983289078

In [110]:
stieltjes(f, 1+0.000001im)


Out[110]:
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 [31]:
x = dual(1,2)
cos(x)


Out[31]:
0.5403023058681398 - 1.682941969615793ɛ

In [32]:
cos(1), -2sin(1)


Out[32]:
(0.5403023058681398, -1.682941969615793)

In [33]:
cos(cos(x))


Out[33]:
0.8575532158463934 + 0.8656973695406469ɛ

In [34]:
cos(cos(1)), 2sin(cos(1))*sin(1)


Out[34]:
(0.8575532158463934, 0.8656973695406469)

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]:
0.16885248872798114 + 0.007397239571120209ɛ

In [36]:
dualpart(ret) - 2*0.0036986197855601063


Out[36]:
-3.469446951953614e-18

The direction can be complex as well:


In [37]:
cos(dual(1.0,im))


Out[37]:
0.5403023058681398 - 0.0im + 0.0ɛ - 0.8414709848078965imɛ

In [38]:
sqrt(dual(-1.0+0im,im))


Out[38]:
0.0 + 1.0im + 0.5ɛ + 0.0imɛ

In [39]:
sqrt(dual(-1.0-0im,-im))


Out[39]:
0.0 - 1.0im + 0.5ɛ + 0.0imɛ

In [40]:
log(dual(0.0,1.0))


Out[40]:
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 [111]:
log(RiemannDual(0.0,1.0))


Out[111]:
(1.0 + 0.0im)log ε + 0.0 + 0.0im

In [42]:
l = log(RiemannDual(0.0,im))


Out[42]:
(1.0 + 0.0im)log ε + 0.0 + 1.5707963267948966im

In [43]:
h = 0.00001
logpart(l)*log(h) + finitepart(l)


Out[43]:
-11.512925464970229 + 1.5707963267948966im

In [44]:
l = stieltjes(f, RiemannDual(1.0,im))

l(h)


Out[44]:
29.593415861168136 - 4.269867111336779im

In [45]:
stieltjes(f, 1.0 + im*h)


Out[45]:
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 [46]:
Γ = Segment(0,1)  Segment(0, im)

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


Out[46]:
-2 -1 0 1 2 -2 -1 0 1 2 y1

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]:
-4.440892098500626e-16 + 0.0im

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]:
-15.133399294132497 + 0.0im

In [50]:
stieltjes(f₂, -0.000001)


Out[50]:
13.575698815971107 + 1.5707948403214076im

But these are cancelled when summed together:


In [51]:
stieltjes(f, -0.000001)


Out[51]:
-1.5577004781613901 + 1.5707948403214076im

This is easy to see with log numbers:


In [52]:
stieltjes(f₁, RiemannDual(0,-1))


Out[52]:
(1.0000000000000009 + 0.0im)log ε + -1.3179021514543992 + 0.0im

In [53]:
stieltjes(f₂, RiemannDual(0,-1))


Out[53]:
(-1.000000000000001 + 0.0im)log ε + -0.23981174200057237 + 1.5707963267948983im

In [54]:
stieltjes(f, RiemannDual(0,-1))


Out[54]:
(-2.220446049250313e-16 + 0.0im)log ε + -1.5577138934549715 + 1.5707963267948983im

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]:
(-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 [56]:
G₁, G₂, G₃, G₄ = first.(pieces(G))

G₄*G₃*G₂*G₁


Out[56]:
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 [57]:
F = G-I
F₁, F₂, F₃, F₄ = first.(pieces(F))

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


Out[57]:
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 [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]:
6.483657258559754e-16

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

In [60]:
norm( F₁*G₂*G₃*G₄ + F₂*G₃*G₄ + F₃*G₄ + F₄)


Out[60]:
6.841415457043595e-18

2017 v 2023 (Future directions.)

President Ivanka Trump.... 😰

Banded operators

Chebyshev singularities lead to banded operators:


In [61]:
H = Hilbert() : ChebyshevWeight()


Out[61]:
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 [62]:
 = Fun(x->exp(-40x^2),ChebyshevWeight())
G = Fun(x->1+exp(-40x^2),Chebyshev())
(I+G)/2 + im/2**H


Out[62]:
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 [63]:
Fun(z -> stieltjes(P(15),z), Legendre(Segment(1+im,2+2im))) |>norm


Out[63]:
8.244707787129415e-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.