In [7]:
using ApproxFun, SingularIntegralEquations, RiemannHilbert, Plots, QuadGK, DualNumbers, ComplexPhasePortrait, LinearAlgebra
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
In [4]:
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[4]:
Associated with the continuous spectrum is the reflection coefficient. This can be calculated:
In [5]:
R = ReflectionCoefficient( x -> 1.2exp(-x^2))
@time R(0.1)
Out[5]:
Behind the scenes, its solving an ODE for $\psi$ on (-∞,0] and $\phi_\pm$ on [0,∞), by truncating the domain:
In [8]:
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\\_-")
Out[8]:
In [9]:
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[9]:
In [10]:
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[10]:
In [11]:
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[11]:
The first step is to expand $R$ into a Chebyshev expansion. We use tFun
which takes advantage of parallelisation:
In [17]:
@time ρ = tFun(R, -5.0..5, 300)
plot(ρ)
Out[17]:
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 [18]:
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 [19]:
@time Φ = transpose(rhsolve(transpose(G), 2*4*200)); # asymptotic to I
Φ = [1 1]*Φ; # asymptotic to [1 1]
It worked!
In [20]:
Φ(0.1+0.0im) ≈ Φ(0.1-0.0im)*G(0.1)
Out[20]:
We can then recover $V$ from Φ
. I'm not going to show that since there's a bug 😳
In [21]:
ψ(0.0)/a - (φ₋(0.0) + ρ(1)*φ₊(0.0))
Out[21]:
In [22]:
Φ(1.0-0.0im)
Out[22]:
In [23]:
Φ(1.0+0.0im)
Out[23]:
In [24]:
R(1)
Out[24]:
In [25]:
ψ(0.0)/a
Out[25]:
In [33]:
(ϕ₋)(0.0)
In [34]:
ϕ₊(0.0)
In [35]:
S = []
Out[35]:
In [36]:
ψ(-10.0)exp(im*1*(-10.0))
Out[36]:
In [37]:
Φ(1-0.0im)
Out[37]:
In [38]:
ϕ₊(0.0)
In [39]:
ϕ₋(1)
In [ ]: