Root-finding

Functions

This is the midpoint RK routine from Chapter 2.

RKODE[G_, x0_,f0_, dx_, Ns_] := Module[{k1,k2,retvals,index,xn,fn},
retvals = Table[{0.0,0.0},{i,1,Ns}];
xn = x0;
fn = f0;
retvals[[1]] = {x0,f0};
For[index = 2, index ≤ Ns, index=index+1,
k1 = dx G[xn,fn];
k2 = dx G[xn + dx/2, fn + k1/2];
fn = fn + k2;
xn = xn + dx;
retvals[[index]] = {xn,fn};
];
Return[retvals];
]

Example for Implementation 3.1

This function takes an energy, Et, solves Schrodinger’s equation for an infinite square well potential using that energy, and returns the value of Psi at the right-hand end-point:

F[Et_] := Module[{Nsteps,G,nvals,PsiE},
Nsteps = 100;
G[x_,f_] := {f[[2]], -f[[3]]  f[[1]],0};
nvals = RKODE[G, 0.0, {0.0,1.0,Et},1.0/(Nsteps-1), Nsteps];
PsiE = nvals[[Length[nvals],2,1]];
Return[PsiE];
]

We can plot the function to see its root structure

G0 = Plot[F[x],{x,0,(5 Pi)^2}, PlotRange->All,PlotStyle->GrayLevel[0.5]]

the roots should be at (n^2)(π^2) for integer n - let’s put circles at the relevant locations (this is Figure 3.5):

rtlist = Table[Point[{ n^2 Pi^2,0}],{n,1,5}];
G1 = Show[Graphics[rtlist]];
Show[G0,G1,AxesOrigin->{0,-.2}]

Bisection - Implementation 3.2

The recursive form of bisection takes a function to bisect, F, and a bracketing [xl, xr] (with F(xl) F(xr) < 0, so there is a root in between) - once |F(xm)| < eps, the function returns - otherwise, it generates a new bracketing of half the original length and calls itself:

Bisect[F_, xl_, xr_, eps_] := Module[{val, vam, xm, retval},
xm = .5(xl + xr);
val = F[xl];
vam = F[xm];
If[Abs[vam] ≤ eps,
Return[xm];
];
If[val vam < 0,
retval = Bisect[F,xl,xm,eps];
,
retval = Bisect[F,xm, xr, eps];
];
Return[retval];
]

We can test it on the function F from above,

E0 = Bisect[F,5.0,12.0,10^(-5)]

This is supposed to be (π^2)

E0 - Pi^2

the accuracy is limited by the numerical solution to the ODE.

Newton’s Method

Newton’s method uses derivative information to estimate the location of a root, and then refines based on that location - it will work if you start “near the root” (whatever that means).

Implementation 3.3

Given F and its derivative DF, Newton’s method starts at x = x0 and updates x according to (3.69) until |F(x)| < eps:

Newton[F_,DF_,x0_,eps_] := Module[{x, deltax,val},
x = x0;
val = F[x];
While[Abs[val] ≥ eps,
deltax = -F[x]/DF[x];
x = x + deltax;
val = F[x];
];
Return[x];
]

We have to provide the function and its derivative

exampleF[x_] := (x-1)(x-2)

exampleDF[x_] := -3 + 2 x

Then we can find the roots, there are two, so we need to be careful to bracket each one separately

Newton[exampleF,exampleDF,.5,10^(-8)]

Newton[exampleF,exampleDF,1.75,10^(-8)]

Implementation 3.4

If we don’t have immediate access to the derivative of F, we can approximate it internally (the “secant” method) - here we must give a delta for the derivative estimate:

nNewton[F_,delta_, x0_, eps_] := Module[{x, deltax, val, DF},
DF[x_] := (F[x + delta] - F[x - delta])/(2. delta);
x = x0;
val = F[x];
While[Abs[val] ≥ eps,
deltax = -F[x]/DF[x];
x = x + deltax;
val = F[x];
];
Return[x];
]

We can find the root of the exampleF without giving the derivative:

nNewton[exampleF,.01,.2,10^(-8)]

or use the new function in the infinite square well case from above, where the derivative of F is simply not available (because F is a numerically defined function):

nNewton[F,.01,5.0,10^(-8)]

Gradients (Figure 3.9)

Here’s an example of a contourplot with an overlaid gradient field - the arrows point towards the maxima.

u[x_,y_] := x^2 - y^2

DU[x_,y_] := {2 x, -2 y}

G0 = VectorPlot[DU[x,y],{x,-1,1},{y,-1,1},VectorStyle->GrayLevel[1]];
G1 = ContourPlot[u[x,y],{x,-1,1},{y,-1,1},Contours->25,ColorFunction->"GrayTones"];
Grout = Show[G1,G0]

Scattering Examples

An implementation of Algorithm 3.1 for two dimensional steepest descent - we must have functions u and du that give back the value of u at (x,y) and the gradient ∇u at (x,y).

Descent[u_, du_, x0_, y0_, eta_,eps_] := Module[{x,y,r},
x = x0;
y = y0;
r= u[x,y];
While[r > eps,
{x, y} = {x,y} - eta du[u,x,y];
r = u[x,y];
];
Return[{x,y}];
]

The function u, for the quantum mechanical scattering cases, will take the real and imaginary parts of the Psi at q = 0 as “x” and “y” - it returns ((Abs[Psip - i Sqrt[Et] Psi])^2) as in (3.64).

Finite Well (Figure 3.10)

Given Et = 2 and a finite well of depth V0 = -10, we can solve numerically and output the relevant combination of Psi (both real and imaginary parts) and its derivative at the right-hand end-point:

uFiniteWell[rPsi0_, iPsi0_ ] := Module[{fval, retval, Psi0,G,Nsteps,Et, Vt},
Et = 2.;
Vt[q_] := -10.;
G[x_,f_] := {f[[2]],-(Et - Vt[x]) f[[1]]};
Nsteps = 1000;
Psi0 = rPsi0 + I iPsi0;
fval = RKODE[G, 0.0,{Psi0, I Sqrt[Et] (2.0 - Psi0)}, 1.0/(Nsteps - 1), Nsteps];
retval = fval[[Length[fval],2,2]] - I Sqrt[Et] fval[[Length[fval],2,1]];
Return[Chop[Abs[retval]^2]];
]

To find the gradient, we use a symmetric difference in each direction with dr set small:

duFiniteWell[u_,rPsi0_, iPsi0_] := Module[{cur, ngrad,dr},
dr = .001;
ngrad = {(u[rPsi0 + dr, iPsi0] - u[rPsi0 - dr,iPsi0])/(2. dr), (u[rPsi0, iPsi0 + dr] - u[rPsi0,iPsi0-dr])/(2. dr)};
Return[ngrad];
]

eps = 10^(-8);
{rPsi0, iPsi0} = Descent[uFiniteWell, duFiniteWell,1.0,2.0,.01, eps];

We can check the value of u given this real and imaginary initial Psi

uFiniteWell[rPsi0, iPsi0]

Now we want to generate the wavefunctions

Psi0 = rPsi0 + I iPsi0;

First we’ll get the numerical portion of the solution:

Et = 2.0;
Vt[x_] := -10.0;
G[x_,f_] := {f[[2]],-(Et - Vt[x]) f[[1]]};
Nsteps = 1000;
fval = RKODE[G, 0.0,{Psi0, I Sqrt[Et] (2.0 - Psi0)}, 1.0/(Nsteps - 1), Nsteps];

Now we can pick up the coefficients we need for the exact, external wavefunctions, the oscillatory pieces:

b = Psi0 -1;
c = fval[[Length[fval],2,1]] Exp[-I Sqrt[Et]];

Psil =  Exp[I Sqrt[Et] q] + b Exp[-I Sqrt[Et] q];
Psir = c Exp[I Sqrt[Et q]];

Finally, plot the probability density for all three regions:

G0 = ListPlot[Table[{fval[[j,1]],Abs[fval[[j,2,1]]]^2},{j,1,Length[fval]}], PlotRange-> {{-2,2},{0,2}}];
G1 = Plot[Abs[Psil]^2,{q,-2,0}, PlotRange->{{-2,2},{0,2}}];
G2 = Plot[Abs[Psir]^2,{q,1,2},PlotRange->{{-2,2},{0,2}}];
Show[{G0,G1,G2}]

And this time, plot the numerical solution using lines, just to avoid the overlapping points:

G0 = ListLinePlot[Table[{fval[[j,1]],Abs[fval[[j,2,1]]]^2},{j,1,Length[fval]}], PlotRange-> {{-2,2},{0,2}},PlotStyle->GrayLevel[0.]];
G1 = Plot[Abs[Psil]^2,{q,-2,0}, PlotRange->{{-2,2},{0,2}},PlotStyle->GrayLevel[0.5]];
G2 = Plot[Abs[Psir]^2,{q,1,2},PlotRange->{{-2,2},{0,2}}, PlotStyle->GrayLevel[0.5]];
Show[{G0,G1,G2}]

Oscillatory Potential (Figure 3.11)

This time, take Et = 2 and an oscillatory V, defined inside u:

uOsc[rPsi0_, iPsi0_ ] := Module[{fval, retval, Psi0,G,Nsteps,Et, Vt},
Et = 2.;
Vt[q_] := 10.0Sin[8 Pi q];
G[x_,f_] := {f[[2]],-(Et - Vt[x]) f[[1]]};
Nsteps = 1000;
Psi0 = rPsi0 + I iPsi0;
fval = RKODE[G, 0.0,{Psi0, I Sqrt[Et] (2.0 - Psi0)}, 1.0/(Nsteps - 1), Nsteps];
retval = fval[[Length[fval],2,2]] - I Sqrt[Et] fval[[Length[fval],2,1]];
Return[Chop[Abs[retval]^2]];
]

duOsc[u_, rPsi0_, iPsi0_] := Module[{cur, ngrad,dr},
dr = .001;
ngrad = {(u[rPsi0 + dr, iPsi0] - u[rPsi0 - dr,iPsi0])/(2. dr), (u[rPsi0, iPsi0 + dr] - u[rPsi0,iPsi0-dr])/(2. dr)};
Return[ngrad];
]

eps = 10^(-8);
{rPsi0, iPsi0} = Descent[uOsc, duOsc,1.0,2.0,.01, eps];

Check the zero:

uOsc[rPsi0, iPsi0]

Generate the initial values and find the numerical solution:

Psi0 = rPsi0 + I iPsi0;
Et = 2.0;
Vt[x_] :=10.0Sin[8 Pi x];
G[x_,f_] := {f[[2]],-(Et - Vt[x]) f[[1]]};
Nsteps = 1000;
fval = RKODE[G, 0.0,{Psi0, I Sqrt[Et] (2.0 - Psi0)}, 1.0/(Nsteps - 1), Nsteps];

Get left and right coefficients:

b = Psi0 -1;
c = fval[[Length[fval],2,1]] Exp[-I Sqrt[Et]];

Psil =  Exp[I Sqrt[Et] q] + b Exp[-I Sqrt[Et] q];
Psir = c Exp[I Sqrt[Et q]];

Probability density for all three regions:

G0 = ListPlot[Table[{fval[[j,1]],Abs[fval[[j,2,1]]]^2},{j,1,Length[fval]}], PlotRange-> {{-2,2},{0,2}}];
G1 = Plot[Abs[Psil]^2,{q,-2,0}, PlotRange->{{-2,2},{0,2}}];
G2 = Plot[Abs[Psir]^2,{q,1,2},PlotRange->{{-2,2},{0,2}}];
Show[{G0,G1,G2}]

And the smoothed version:

G0 = ListLinePlot[Table[{fval[[j,1]],Abs[fval[[j,2,1]]]^2},{j,1,Length[fval]}], PlotRange-> {{-2,2},{0,2}},PlotStyle->GrayLevel[0.]];
G1 = Plot[Abs[Psil]^2,{q,-2,0}, PlotRange->{{-2,2},{0,2}},PlotStyle->GrayLevel[0.5]];
G2 = Plot[Abs[Psir]^2,{q,1,2},PlotRange->{{-2,2},{0,2}}, PlotStyle->GrayLevel[0.5]];
Show[{G0,G1,G2}]