Chapter 10
The Method of Inverse Scattering

10.2 The Inverse Scattering Transform

Example 10.1 (Figure 10.2)

Mathematica can be used to examine the analytic structure of the transmission coefficient in a simple example that is well-known to students of introductory quantum mechanics courses, namely, the square-well potential: u(x)=0 if |x|>a and u(x)=-V if |x|<a. The potential well has width a, depth V, and is centered at the
origin.  The calculation of the transmission coefficient for this potential is discussed, for example, by Landau and Lifshitz (1965).

In Figure 10.2 the transmission coefficient a[k] is shown in the upper half of the complex k plane as both surface and contour plots, with the well width and well depth chosen to support two bound states.  The singularities in the transmission coefficient for points on the imaginary k axis are evident in both the surface and the contour representation.  These singularities occur at values of ik that correspond to the bound states of the potential.

The two panels of Figure 10.2 are obtained by constructing the transmission coefficient and then displaying the behavior of this quantity as a function of the complex wave vector using Plot3D and ContourPlot.

a[E_,V_]=1/(1+((V Sin[Sqrt[E-V]])^2/(4E(E-V))))

Plot3D[Abs[a[(Rek+I Imk)^2,-(3Pi/2)^2]],{Rek,-1,1},
     {Imk,0,5},PlotPoints->70,PlotRange->{0,5},
     BoxRatios->{1,1,1},ViewPoint->{-2.4,-1.3,2},
     Shading->False,AxesLabel->{Re,Im,a},
     ClipFill->None]

ContourPlot[Abs[a[(Rek+I Imk)^2,-(3Pi/2)^2]],
     {Rek,-1,1},{Imk,0,5},PlotPoints->80,
     Contours->20,ContourShading->False,
     Ticks->{Range[-1,1,0.5],Range[0,5,1]}]

10.3 Solution of the Scattering Problem

The computer program shown below is taken from the book by Richard Crandall (Mathematica for the Sciences, Addison-Wesley, 1991).  It is robust enough to handle a variety of smooth initial conditions, but is not stable enough to handle a step function.  The parameters are:
     
     kmax, the number of spatial grid points
     kplot, the grid points used for plotting the solution
     x0, the "center" of the initial function
     dt, the time step of the integration
     dx, the spatial increment of the grid points
     iter, the number of iterations, i.e. iter*dt is the time span of the integration
     view, the frequency with which the solution is plotted, i.e. every viewth
              iteration
     F[x_], the initial condition of the Korteweg-de Vries equation

Figure 10.4

F[x_] = 2(Sech[x])^2;
kmax  = 200;
kplot = 175;
x0    = 150;
dt    = 0.002;
dx    = 0.2;
iter  = 225;
view  = 75;
uPast = Table[F[(k-x0)dx],{k,1,kmax}];
uPres = uPast;
uFutu = uPres;
Do[uPres[[j]] = uPast[[j]] +
    dt (6 uPast[[j]](uPast[[j+1]]-uPast[[j-1]])/(2 (dx)) -
    (uPast[[j+2]] - 2 uPast[[j+1]] + 2 uPast[[j-1]] -
    uPast[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
Do[Do[uFutu[[j]] = uPast[[j]] +
    2 dt (6 uPres[[j]](uPres[[j+1]]-uPres[[j-1]])/(2 (dx)) -
    (uPres[[j+2]] - 2 uPres[[j+1]] + 2 uPres[[j-1]] -
    uPres[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
uPast = uPres;
uPres = uFutu;
uPlot = Table[uFutu[[j]],{j,1,kplot}];
If[Mod[n,view]==0,
      ListPlot[uPlot, PlotRange->{-1,2},
      AxesLabel->{x,u},Ticks->{Range[0,150,50],
      Range[-1,2,1]},PlotJoined->True]
      ],{n,0,iter}]

Figure 10.5

F[x_] = -1(Sech[x])^2;
kmax  = 200;
kplot = 175;
x0    = 100;
dt    = 0.002;
dx    = 0.2;
iter  = 750;
view  = 250;
uPast = Table[F[(k-x0)dx],{k,1,kmax}];
uPres = uPast;
uFutu = uPres;
Do[uPres[[j]] = uPast[[j]] +
    dt (6 uPast[[j]](uPast[[j+1]]-uPast[[j-1]])/(2 (dx)) -
    (uPast[[j+2]] - 2 uPast[[j+1]] + 2 uPast[[j-1]] -
    uPast[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
Do[Do[uFutu[[j]] = uPast[[j]] +
    2 dt (6 uPres[[j]](uPres[[j+1]]-uPres[[j-1]])/(2 (dx)) -
    (uPres[[j+2]] - 2 uPres[[j+1]] + 2 uPres[[j-1]] -
    uPres[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
uPast = uPres;
uPres = uFutu;
uPlot = Table[uFutu[[j]],{j,1,kplot}];
If[Mod[n,view]==0,
      ListPlot[-uPlot, PlotRange->{-0.2,1},
      AxesLabel->{x,u},Ticks->{Range[0,150,50],
      Range[-0.2,1,0.2]},PlotJoined->True]
      ],{n,0,iter}]

Figure 10.6

F[x_] = -2(Sech[x])^2;
kmax  = 200;
kplot = 175;
x0    = 100;
dt    = 0.002;
dx    = 0.2;
iter  = 750;
view  = 250;
uPast = Table[F[(k-x0)dx],{k,1,kmax}];
uPres = uPast;
uFutu = uPres;
Do[uPres[[j]] = uPast[[j]] +
    dt (6 uPast[[j]](uPast[[j+1]]-uPast[[j-1]])/(2 (dx)) -
    (uPast[[j+2]] - 2 uPast[[j+1]] + 2 uPast[[j-1]] -
    uPast[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
Do[Do[uFutu[[j]] = uPast[[j]] +
    2 dt (6 uPres[[j]](uPres[[j+1]]-uPres[[j-1]])/(2 (dx)) -
    (uPres[[j+2]] - 2 uPres[[j+1]] + 2 uPres[[j-1]] -
    uPres[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
uPast = uPres;
uPres = uFutu;
uPlot = Table[uFutu[[j]],{j,1,kplot}];
If[Mod[n,view]==0,
      ListPlot[-uPlot, PlotRange->{0,2},
      AxesLabel->{x,u},Ticks->{Range[0,150,50],
      Range[0,2,1]},PlotJoined->True]
      ],{n,0,iter}]

Figure 10.7

F[x_] = -4(Sech[x])^2;
kmax  = 200;
kplot = 185;
x0    = 100;
dt    = 0.002;
dx    = 0.2;
iter  = 750;
view  = 250;
uPast = Table[F[(k-x0)dx],{k,1,kmax}];
uPres = uPast;
uFutu = uPres;
Do[uPres[[j]] = uPast[[j]] +
    dt (6 uPast[[j]](uPast[[j+1]]-uPast[[j-1]])/(2 (dx)) -
    (uPast[[j+2]] - 2 uPast[[j+1]] + 2 uPast[[j-1]] -
    uPast[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
Do[Do[uFutu[[j]] = uPast[[j]] +
    2 dt (6 uPres[[j]](uPres[[j+1]]-uPres[[j-1]])/(2 (dx)) -
    (uPres[[j+2]] - 2 uPres[[j+1]] + 2 uPres[[j-1]] -
    uPres[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
uPast = uPres;
uPres = uFutu;
uPlot = Table[uFutu[[j]],{j,1,kplot}];
If[Mod[n,view]==0,
      ListPlot[-uPlot, PlotRange->{-0.5,5},
      AxesLabel->{x,u},Ticks->{Range[0,150,50],
      Range[0,5,1]},PlotJoined->True]
      ],{n,0,iter}]

Figure 10.8

F[x_] = -6(Sech[x])^2;
kmax  = 250;
kplot = 175;
x0    = 100;
dt    = 0.001;
dx    = 0.2;
iter  = 750;
view  = 250;
uPast = Table[F[(k-x0)dx],{k,1,kmax}];
uPres = uPast;
uFutu = uPres;
Do[uPres[[j]] = uPast[[j]] +
    dt (6 uPast[[j]](uPast[[j+1]]-uPast[[j-1]])/(2 (dx)) -
    (uPast[[j+2]] - 2 uPast[[j+1]] + 2 uPast[[j-1]] -
    uPast[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
Do[Do[uFutu[[j]] = uPast[[j]] +
    2 dt (6 uPres[[j]](uPres[[j+1]]-uPres[[j-1]])/(2 (dx)) -
    (uPres[[j+2]] - 2 uPres[[j+1]] + 2 uPres[[j-1]] -
    uPres[[j-2]])/(2 (dx)^3)),{j,3,kmax-2}];
uPast = uPres;
uPres = uFutu;
uPlot = Table[uFutu[[j]],{j,1,kplot}];
If[Mod[n,view]==0,
      ListPlot[-uPlot, PlotRange->{0,10},
      AxesLabel->{x,u},Ticks->{Range[0,150,50],
      Range[0,10,2]},PlotJoined->True]
      ],{n,0,iter}]

10.4 Lax's Method

To carry out the algebraic steps in (10.115) with Mathematica, the operators  L and B are first defined:

f:=Psi[x,y]

g:=u[x,y]

LaxL[f_]:=-D[f,{x,2}] + g f

LaxB[f_]:=-4 D[f,{x,3}] + 6 g D[f,x] + 3 f D[g,x]

The operations in equations (10.115) are now carried out with the command

Simplify[LaxB[LaxL[f]]-LaxL[LaxB[f]]-
     Dt[LaxL[f],y,Constants->{D[f,{x,2}],f,x}]]

which, when set equal to zero, yields (10.119).  Notice that in calculating the
left-hand side of this equation, the operation Dt was applied to L[Psi] and only u was considered as being t-dependent.