Numerical Inverse Scattering Transform for Solving the Nonlinear - - PowerPoint PPT Presentation

Numerical Inverse Scattering Transform for Solving the Nonlinear Schr¨

• dinger Equation
• A. Aric`
• ,
• G. Rodriguez,
• S. Seatzu

Department of Mathematics and Computer Science University of Cagliari · Italy Conference on Applied Inverse Problems July 20-24, 2009 · Vienna · Austria

Outline

Nonlinear Schr¨

• dinger equation and Inverse Scattering Transform

Nonlinear Schr¨

• dinger equation (NLS)

Numerical method Issues related to the application of the IST to the NLS Numerical solution to the inverse scattering problem Time evolution of the scattering data Solution to the Marchenko systems Numerical experiments One soliton Truncated one soliton

Numerical solution to the cubic NLS - AIP09 · Vienna slide 2/10

Nonlinear Schr¨

• dinger equation (NLS)

We are interested in the IVP for the cubic Schr¨

• dinger equation (NLS)

   iqt = qxx + 2q|q|2, x ∈ R, t > 0, q(x; 0) given x ∈ R, where q = q(x; t); we follow the path of the Inverse Scattering Transform (IST): initial potential q(x; 0)

NLS: iqt = qxx + 2q|q|2

• direct

scattering

• initial scattering data

S(α; 0)

time evolution

• f scattering data:

iSt = ±4Sαα

• evolved potential

q(x; t)

• inverse

scattering

evolved scattering data S(α; t)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 3/10

Nonlinear Schr¨

• dinger equation (NLS)

We are interested in the IVP for the cubic Schr¨

• dinger equation (NLS)

   iqt = qxx + 2q|q|2, x ∈ R, t > 0, q(x; 0) given x ∈ R, where q = q(x; t); we follow the path of the Inverse Scattering Transform (IST): initial potential q(x; 0)

NLS: iqt = qxx + 2q|q|2

• direct

scattering

• initial scattering data

S(α; 0)

time evolution

• f scattering data:

iSt = ±4Sαα

• evolved potential

q(x; t)

• inverse

scattering

evolved scattering data S(α; t)

“The IST solves exactly the IVP for the NLS; all the physical properties are then preserved.”

Numerical solution to the cubic NLS - AIP09 · Vienna slide 3/10

Numerical method

Here we validate the last two step of the procedure. A further work will be the implementation of the direct scattering, thanks to a result due to van der Mee.

Proposed algorithm

• 1. assume that we know S(λ; 0); choose {t1, . . . , tk}
• 2. for each ti:

◮ choose a discretization step h and sample S(λ; 0) with step h: Sh(0) ◮ apply a numerical method to Sh(0) and compute Sh(ti) ◮ apply a numerical method to Sh(ti) and compute [twice] q(0; ti) ◮ check if q(0; ti) is accurate:

YES compute q(xk; ti) for all xk NO halve h and restart the computation for the same ti

Numerical solution to the cubic NLS - AIP09 · Vienna slide 4/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

Scattering data

◮ “reflection coefficients” from the left and from the right: L(λ), R(λ); ◮ “transmission coefficients”: T(λ); ◮ “bound states” λj ∈ C+, and their “norming constants” Γlj, Γrj.

In the numerical experiments for the Inverse Scattering, we assume to know the initial scattering data and use them to compute the potential.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

Time evolution of the relevant scattering data

L(λ)

e−4iλ2t

• R(λ)

e+4iλ2t

• λj

1

• Γlj

e+4iλ2t

• Γrj

e−4iλ2t

• L(λ; t)

R(λ; t) λj(t) Γlj(t) Γrj(t) The scattering data determines the Marchenko integral kernels, i.e. the kernels of the integral system related to the inverse scattering step.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

Kernels of the Marchenko integral systems

Marchenko kernel from the left Marchenko kernel from the right Ωl(α) = ˆ R(α) +

• j

Γljeiλjα

i(Ωl)t = 4(Ωl)αα

• Ωr(α) = ˆ

L(α) +

• j

Γrjeiλjα

i(Ωr)t = −4(Ωr)αα

• Ωr(α; t) = ˆ

R(α; t) +

• j

Γrj(t)eiλjα Ωr(α; t) = ˆ L(α; t) +

• j

Γrj(t)eiλjα

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

Marchenko integral systems:

For a fixed value t 0, and for any value of x ∈ R, solve: (B(R, R+; t)) left        Bl1(x, γ; t) − +∞ Ωl(2x + γ + β; t)Bl2(x, β; t) dβ = 0 +∞ Ωl(2x + γ + β; t)Bl1(x, β; t) dβ + Bl2(x, γ; t) = −Ωl(2x + γ; t) and/or right        Br1(x, γ; t) + +∞ Ωr(−2x + γ + β; t)Br2(x, β; t) dβ = 0, − +∞ Ωr(−2x + γ + β; t)Br1(x, β; t) dβ + Br2(x, γ; t) = Ωr(−2x + γ; t),

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Issues related to the application of the IST to the NLS

◮ eigenvalue system + FT + nonlinear least square problem; ◮ linear evolution problem; ◮ two systems of two Marchenko integral equations

The Marchenko system returns the solution:

Potential q(x; t) =    2Bl2(x, 0+; t), x 0, −2Br2(x, 0+; t), x 0; Density energy:          Bl1(x, 0+; t) = −1 2 +∞

x

|q(y; t)|2 dy, Br1(x, 0+; t) = −1 2 x

−∞

|q(y; t)|2 dy,

Numerical solution to the cubic NLS - AIP09 · Vienna slide 5/10

Time evolution of the scattering data

◮ The time evolution of each bound state term is straightforward:

Γljeiλjα →

• Γljeiλ2

j ti

eiλjα, Γljeiλjα →

• Γljeiλ2

j ti

eiλjα.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 6/10

Time evolution of the scattering data

◮ The time evolution of each bound state term is straightforward:

Γljeiλjα →

• Γljeiλ2

j ti

eiλjα, Γljeiλjα →

• Γljeiλ2

j ti

eiλjα.

◮ About the time evolution of the reflection coefficients L(λ), R(λ):

L(λ) e−4iλ2ti

• L(λ; ti)

F ˆ L(α; ti) R(λ) e4iλ2ti

• R(λ; ti)

F ˆ R(α; ti)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 6/10

Time evolution of the scattering data

◮ The time evolution of each bound state term is straightforward:

Γljeiλjα →

• Γljeiλ2

j ti

eiλjα, Γljeiλjα →

• Γljeiλ2

j ti

eiλjα.

◮ About the time evolution of the reflection coefficients L(λ), R(λ):

L(λ) e−4iλ2ti

• F−1

ˆ L(α)

• iˆ

Lt = −4ˆ Lαα

• L(λ; ti)

F ˆ L(α; ti) R(λ) e4iλ2ti

• F−1

ˆ R(α)

• iˆ

Rt = 4ˆ Rαα

• R(λ; ti)

F ˆ R(α; ti)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 6/10

Time evolution of the scattering data

◮ The time evolution of each bound state term is straightforward:

Γljeiλjα →

• Γljeiλ2

j ti

eiλjα, Γljeiλjα →

• Γljeiλ2

j ti

eiλjα.

◮ About the time evolution of the reflection coefficients L(λ), R(λ):

L(λ) e−4iλ2ti

• F−1

ˆ L(α)

• iˆ

Lt = −4ˆ Lαα

• L(λ; ti)

F ˆ L(α; ti) R(λ) e4iλ2ti

• F−1

ˆ R(α)

• iˆ

Rt = 4ˆ Rαα

• R(λ; ti)

F ˆ R(α; ti)

• 1. [α0, αN] support of ˆ

R(α), i.e. at time t = 0;

• 2. approximate ˆ

R(α; ti) on [α0, αN]: apply IFFT; integrate (multiply by e4iλ2ti); apply FFT;

• 3. check if ˆ

R(α; ti) is negligible at the boundary; (4.) enlarge [α0, αN] and recompute ˆ R(α; ti) with the same stepsize.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 6/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

left        Bl1(x, γ; t) − +∞ Ωl(2x + γ + β; t)Bl2(x, β; t) dβ = 0 +∞ Ωl(2x + γ + β; t)Bl1(x, β; t) dβ + Bl2(x, γ; t) = −Ωl(2x + γ; t) right        Br1(x, γ; t) + +∞ Ωr(−2x + γ + β; t) Br2(x, β; t) dβ = 0, − +∞ Ωr(−2x + γ + β; t)Br1(x, β; t) dβ + Br2(x, γ; t) = Ωr(−2x + γ; t), 1) conjugate both equations; 2) change of sign in the first integral and in the second equation; 3) replace x by −x

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

left        Bl1(x, γ; t) − +∞ Ωl(2x + γ + β; t)Bl2(x, β; t) dβ = 0 +∞ Ωl(2x + γ + β; t)Bl1(x, β; t) dβ + Bl2(x, γ; t) = −Ωl(2x + γ; t) right        Br1(−x, γ; t) − +∞ Ωr(2x + γ + β; t) (−Br2(−x, β; t)) dβ = 0, +∞ Ωr(2x + γ + β; t)Br1(−x, β; t) dβ + (−Br2(−x, γ; t)) = −Ωr(2x + γ; t), 1) conjugate both equations; 2) change of sign in the first integral and in the second equation; 3) replace x by −x

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

left        Bl1(x, γ; t) − +∞ Ωl(2x + γ + β; t)Bl2(x, β; t) dβ = 0 +∞ Ωl(2x + γ + β; t)Bl1(x, β; t) dβ + Bl2(x, γ; t) = −Ωl(2x + γ; t) Given xk and t, we collocate the left Marchenko system in each point of the grid {βi = ih, i = 0, . . . , N}, and approximate the integrals by appliyng the repeated Simpson’s quadrature rule in the same points.

◮ h is the stepsize used in the time evolution routine ◮ N is the number of points where ˆ

R(α; ti) has been approximated

◮ xk belongs to the same interval [ α0 2 , αN 2 ], with half stepsize

xk = αk 2 = α0 2 + h 2k ⇒ 2xk + γi + βj = αk+i+j

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

left        Bl1(x, γ; t) − +∞ Ωl(2x + γ + β; t)Bl2(x, β; t) dβ = 0 +∞ Ωl(2x + γ + β; t)Bl1(x, β; t) dβ + Bl2(x, γ; t) = −Ωl(2x + γ; t) Given xk and t, we collocate the left Marchenko system in each point of the grid {βi = ih, i = 0, . . . , N}, and approximate the integrals by appliyng the repeated Simpson’s quadrature rule in the same points. left              Bl1(xk, βi; t) −

N

• j=0

Ωl(αk+i+j; t)djBl2(xk, βj; t) = 0

N

• j=0

Ωl(αk+i+j; t)djBl1(xk, βi; t) + Bl2(xk, βi; t) = −Ω(αk+i; t)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

The support of Ωl(α; t) is essentially (−∞, αN], so we put Ωl(α; t) = 0 whenever α > αt

• N. In matrix notation, the previous sistem takes the form

I −H∗

l D

HlD I bl1 bl2

• =

−ω

• where, for i, j = 0, . . . , N,

◮ D = h 3 diag(1, 4, 2, 4, 2, . . . ) ◮ (Hl)ij = Ωl(αk+i+j; t)

(upper-left triangular + low rank)

◮ (bl1)i = B1l(xk, βi; t) ◮ (bl2)i = B2l(xk, βi; t)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

Solution to the Marchenko systems

◮ one routine to solve both the Marchenko systems ◮ discretization of the “left” system by the Nystr¨

• m method

◮ numerical method to solve the linear system

If the reflection coefficient vanishes (L(λ) = R(λ) = 0), the system matrix is “Identity + Low Rank”. Low = twice the number of bound states.

◮ we have an exact formula for the inverse ◮ we can compute the limit of the solution, wrt h, and get the

analytical solution of the Marchenko systems

◮ this has been implemented, and improves the overall results

(precision, time complexity) If the reflection coefficients do not vanish, we apply the CG to the Schur system, i.e. (I + HDH∗D)bl2 = −ω, but in a symmetryzed version: (D + DHDH∗D)bl2 = −Dω. Since the B(x, γ; t) are (. . . ) regular, we can also reduce the number of required iterations.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 7/10

One soliton

◮ parameters:

η, ξ, x0, φ ∈ R, η > 0

◮ initial potential:

q(x; 0) = −2iηe−i(2ξx+φ) sech(x0 − 2ηx)

◮ scattering data:

λ1 = iη − ξ, Γl1 = 2iηex0−iφ Γr1 = −2iηe−x0+iφ T(λ) = λ + λ1 λ − λ1 , L(λ) = 0 R(λ) = 0,

◮ solution to the NLS:

q(x; t) = −2iηe−i(2ξx+4(η2−ξ2)t+φ) sech(x0 − 2ηx + 8ηξt) The amplitude is 2η and the peak point at time t is x0

2η + 4ξt.

Numerical solution to the cubic NLS - AIP09 · Vienna slide 8/10

One soliton

Parameters: x ∈ [−12, 7], η ∈ {1, 3, 6}, ξ = − 1

3,

x0 = 1, φ = π

4 .

IST (h ≃ 0.075) FSS (∆t ≃ 2.6 · 10−4)

−10−5 0 5 2 4 0.5 1 1.5 |u| −10−5 0 5 2 4 5 x 10

−15

|e| −10−5 0 5 2 4 −1 1 Re u −10−5 0 5 2 4 −1 1 Im u −10−5 0 5 2 4 0.5 1 1.5 |u| −10−5 0 5 2 4 2 4 x 10

−6

|u−utrue| −10−5 0 5 2 4 −1 1 Re u −10−5 0 5 2 4 −1 1 Im u

Numerical solution to the cubic NLS - AIP09 · Vienna slide 8/10

One soliton

Parameters: x ∈ [−12, 7], η ∈ {1, 3, 6}, ξ = − 1

3,

x0 = 1, φ = π

4 .

IST (h ≃ 0.075) FSS (∆t ≃ 2.6 · 10−4)

−10−5 0 5 2 4 2 4 6 |u| −10−5 0 5 2 4 0.5 1 x 10

−13

|e| −10−5 0 5 2 4 −5 5 Re u −10−5 0 5 2 4 −5 5 Im u −10−5 0 5 2 4 2 4 6 |u| −10−5 0 5 2 4 1 x 10

−4

|u−utrue| −10−5 0 5 2 4 −5 5 Re u −10−5 0 5 2 4 −5 5 Im u

Numerical solution to the cubic NLS - AIP09 · Vienna slide 8/10

One soliton

Parameters: x ∈ [−12, 7], η ∈ {1, 3, 6}, ξ = − 1

3,

x0 = 1, φ = π

4 .

IST (h ≃ 0.075) FSS (∆t ≃ 2.6 · 10−4)

−10−5 0 5 2 4 2 4 6 8 10 |u| −10−5 0 5 2 4 0.5 1 x 10

−12

|e| −10−5 0 5 2 4 −10 10 Re u −10−5 0 5 2 4 −10 10 Im u −10−5 0 5 2 4 2 4 6 8 10 |u| −10−5 0 5 2 4 0.2 |u−utrue| −10−5 0 5 2 4 −10 10 Re u −10−5 0 5 2 4 −10 10 Im u

Numerical solution to the cubic NLS - AIP09 · Vienna slide 8/10

Truncated one-soliton

We introduce a discontinuity in µ: qµ(x) =

• 0,

x < µ q(x; 0), x µ (1)

◮ all the scattering data can be computed analitically ◮ the reflection coefficients do not vanish ◮ they are discontinuous ◮ parameters: η = 1,

ξ = 0, x0 = − ln(2), φ = π

2 ◮ there is one bound state if µ − 1 2 ln(2); none if µ > − 1 2 ln(2) ◮ numerical experiments with µ ∈ {−1, +1}

Numerical solution to the cubic NLS - AIP09 · Vienna slide 9/10

Truncated one-soliton

We introduce a discontinuity in µ=-1 qµ(x) =

• 0,

x < µ q(x; 0), x µ (1) IST & FSS (h ≃ 0.075) IST (t = 0)

−20 −15 −10 −5 5 10 0.5 1 1.5 |q(x;1)| fss ist −20 −10 10 0.5 1 1.5 Re q(x;1) −20 −10 10 −0.1 0.1 0.2 Im q(x;1) −20 −15 −10 −5 5 10 0.5 1 1.5 |q(x;0)| ist true −20 −15 −10 −5 5 10 0.005 0.01 0.015 |q(x;0)−qtrue(x;0)| Numerical solution to the cubic NLS - AIP09 · Vienna slide 9/10

Truncated one-soliton

We introduce a discontinuity in µ=+1 qµ(x) =

• 0,

x < µ q(x; 0), x µ (1) IST & FSS (h ≃ 0.075) IST (t = 0)

−20 −15 −10 −5 5 10 0.01 0.02 0.03 0.04 |q(x;1)| fss ist −20 −10 10 −0.04 −0.02 0.02 Re q(x;1) −20 −10 10 −0.02 0.02 0.04 Im q(x;1) −20 −15 −10 −5 5 10 0.1 0.2 |q(x;0)| ist true −20 −15 −10 −5 5 10 1 2 x 10

−4

|q(x;0)−qtrue(x;0)| Numerical solution to the cubic NLS - AIP09 · Vienna slide 9/10

References

Gardner, Greene, Kruskal, Miura Method for Solving the Korteweg-deVries Equation

• Phys. Rev. Lett. 19-19 (1967), pp. 1095–1097

Zakharov, Shabat Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media ˘

• Z. `
• Eksper. Teoret. Fiz. 61-1 (1971), pp 118–134

Taha, Ablowitz Analytical and numerical aspects of certain nonlinear evolution

• equations. II. Numerical, nonlinear Schr¨
• dinger equation
• J. Comput. Phys. 55-2 (1984), pp 203–230

van der Mee Direct and inverse scattering for skewselfadjoint Hamiltonian systems in Current trends in operator theory and its applications (2004) Ablowitz, Prinari, Trubatch Discrete and continuous nonlinear Schr¨

• dinger systems

Cambridge University Press (2004)

Numerical solution to the cubic NLS - AIP09 · Vienna slide 10/10