In [27]:

using ApproxFun, SingularIntegralEquations, RiemannHilbert, Plots, QuadGK, DualNumbers, ComplexPhasePortrait
import ApproxFunBase: SequenceSpace, BasisFunctional, ℓ⁰, SpaceOperator, piece, pieces, npieces
import ApproxFun: eigs
import ApproxFunOrthogonalPolynomials: Recurrence
import RiemannHilbert: RiemannDual, logpart, finitepart, istieltjes, LogNumber
import SingularIntegralEquations: ⁺, ⁻
using RiemannHilbert.KdV



## The direct scattering transform

We first consider calculation of discrete spectrum of Schrödinger: $$\psi'' + V(x) \psi = \lambda \psi$$



In [196]:

V = x -> 1.2exp(-x^2)
d = -10..10 # computational domain
D = Derivative()

λ,Q = eigs(Dirichlet(d), D^2 + Fun(V,d), 500)
scatter(λ, zero.(λ); xlims=(-2,2))




Out[196]:

-2

-1

0

1

2

0.00

0.25

0.50

0.75

1.00

y1



Associated with the continuous spectrum is the reflection coefficient. This can be calculated:



In [218]:

R = ReflectionCoefficient(V)
R(0.1)




Out[218]:

-0.9651837057640016 + 0.06801873361190361im



Behind the scenes, its solving an ODE for $\psi$ on (-∞,0] and $\phi_\pm$ on [0,∞), by truncating the domain:



In [214]:

M = 10
V = x -> 1.2exp(-x^2)
V₋,V₊ = Fun(V, (-M)..0), Fun(V, 0..M)

k = 1
@time ψ = [ivp(); D^2  + (V₋ + k^2)] \ [exp(-im*k*(-M)), -im*k*(exp(-im*k*(-M))), 0.0]

F = qr([rdirichlet(space(V₊)); rneumann(); D^2  + (V₊ + k^2)])
φ₊ = F \ [exp(im*k*M), im*k*(exp(im*k*M)), 0.0]
φ₋ = F \ [exp(-im*k*M), -im*k*(exp(-im*k*M)), 0.0]

plot(imag(ψ); legend=:bottomleft, label="\\psi", title="k = $k, imaginary part") plot!(imag(φ₊); label="\\phi\\_+") plot!(imag(φ₋); label="\\phi\\_-")   0.015890 seconds (161.25 k allocations: 6.126 MiB) Out[214]: -10 -5 0 5 10 -1.0 -0.5 0.0 0.5 1.0 k = 1, imaginary part ψ φ _ + φ _ -   In [215]: a,b = [φ₋(0) φ₊(0); φ₋'(0) φ₊'(0) ] \ [ψ(0); ψ'(0)] plot(imag(ψ); legend=:bottomleft, label="\\psi", title="k =$k, imaginary part")
plot!(imag(a*φ₋ + b*φ₊); label="a\\phi\\_+ + b\\phi\\_-")




Out[215]:

-10

-5

0

5

10

-1.0

-0.5

0.0

0.5

1.0

k = 1, imaginary part

ψ

a

φ

_

+ b

φ

_

-

+




In [220]:

q = Q[findmax(λ)[2]]
q = q/norm(q)
plot(x -> -V(x), -10,10; label="potential",  legend=:bottomright, ylims=(-2.5,2.5), linestyle=:dot)
plot!(-findmax(λ)[1]+abs2(q); label="eigenfunction")
plot!(abs2(ψ+a*φ₋ + b*φ₊)+k; label="wave")




Out[220]:

-10

-5

0

5

10

-2

-1

0

1

2

potential

eigenfunction

wave



## The inverse scattering transform

With $R$ we have han efficient way of calculating reflection coefficients. We now see how to solve the Riemann–Hilbert problem, but lets simplify things so that we only have continuous spectrum:



In [221]:

V = x -> 0.1exp(-x^2)
d = -10..10 # computational domain
D = Derivative()

λ,Q = eigs(Dirichlet(d), D^2 + Fun(V,d), 500)
scatter(λ, zero.(λ); xlims=(-2,2))




Out[221]:

-2

-1

0

1

2

0.00

0.25

0.50

0.75

1.00

y1



The first step is to expand $R$ into a Chebyshev expansion. We use tFun which takes advantage of parallelisation:



In [222]:

@time ρ = tFun(R, -5.0..5, 300)
plot(ρ)




5.145820 seconds (73.69 M allocations: 2.392 GiB, 13.47% gc time)

Out[222]:

-4

-2

0

2

4

-1.00

-0.75

-0.50

-0.25

0.00

y1

y2



For any $t$ and $x$ we can construct the jump function, but we would need to deform for large $t$ and $x$. Lets keep life simple and take $t = x = 0$, in which case we get the following jump:



In [223]:

k = Fun(identity, space(ρ))
G = [1-abs2.(ρ) -conj.(ρ);
ρ           1.0];



We can now find the solution to $\Phi_+ = \Phi_- G$. RiemannHilbert.jl uses the transpose version ($\Phi_+ = G \Phi_-$) so we transpose twice:



In [224]:

@time Φ = transpose(rhsolve(transpose(G), 2*4*200)); # asymptotic to I
Φ = [1 1]*Φ; # asymptotic to [1 1]




0.615777 seconds (2.02 M allocations: 282.682 MiB)



It worked!



In [225]:

Φ(0.1+0.0im) ≈ Φ(0.1-0.0im)*G(0.1)




Out[225]:

true



We can then recover $V$ from Φ. I'm not going to show that since there's a bug 😳



In [234]:

ψ(0.0)/a - (φ₋(0.0) + ρ(1)*φ₊(0.0))




Out[234]:

-1.1102230246251565e-16 + 1.1657341758564144e-15im




In [227]:

Φ(1.0-0.0im)




Out[227]:

1×2 Array{Complex{Float64},2}:
1.05043-0.372028im  0.942763+0.332064im




In [230]:

Φ(1.0+0.0im)




Out[230]:

1×2 Array{Complex{Float64},2}:
0.942763-0.332064im  1.05043+0.372028im




In [219]:

R(1)




Out[219]:

-0.07910491062180808 + 0.06606135311851252im




In [232]:

ψ(0.0)/a




Out[232]:

0.6642329847792622 + 0.3535499072338503im




In [231]:

(ϕ₋)(0.0)




Out[231]:

0.9732652639976271 + 0.027633169815842518im




In [233]:

ϕ₊(0.0)




Out[233]:

0.9732652639976271 - 0.027633169815842518im




In [ ]:

S = []




In [164]:

ψ(-10.0)exp(im*1*(-10.0))




Out[164]:

1.0000000000000007 + 5.551115123125783e-17im




In [146]:

Φ(1-0.0im)




Out[146]:

2×2 Array{Complex{Float64},2}:
1.03402-0.113411im  -0.0808899+0.20814im
0.0090288-0.286308im     1.00923+0.134907im




In [158]:

ϕ₊(0.0)




Out[158]:

0.7008654655943309 - 0.28472652028393397im




In [151]:

ϕ₋(1)




Out[151]:

0.5383978200907342 - 0.7947809510658824im




In [ ]: