Matrix Numerov Function which can be Generalized to Many Potentials

Ideas
, , ,
1

I have written a function that can implement the Matrix Numerov method in a similar manner to Pillai, et al at UW-Madison, except that this script can quickly be generalized to other potential functions and run in two lines of code. I have tested it on the V potential which can also be solved by Airy functions. Thus, you can accurately approximate eigenvalues and eigenvectors for a variety of potentials.

First, a Derivation:

The Numerov algorithm is for solving differential equations of the form

y''(x)=-g(x)*y(x) + s(x)

Invoke a Taylor expansion of y(x) about x_0:

y(x)= y(x_0) +(x-x_0)\; y'(x_0) +\frac{(x-x_0)^2 }{2!} y''(x_0) +\frac{(x-x_0)^3}{3!} y'''(x_0) +\frac{(x-x_0)^4}{4!} y^{iv}(x_0)+...

Rewrite, letting h equal the distance from point x_0 to x such that y(x)=y(x0+h):

y(x_0+h) =y(x_0) +h\;y'(x_0) +\frac{h^2}{2!} y''(x_0) +\frac{h^3}{3!}\;y'''(x_0) +\frac{h^4}{4!} y^{iv}(x_0)+...

Define a grid of x values and let (x0+h)=x_{n+1}
and x0=x_n.
Then a step forward (toward positive x) is:

y_{n+1} =y_{n}+ h\;y'(x_{n} ) +\frac{h^2}{2!}\;y''(x_{n}) +\frac{h^3}{3!}\;y'''(x_{n})+ \frac{h^4}{4!} y^{iv}(x_n)+...

That was a step forwards. For a step backwards we replace h with -h:

y_{n-1} =y_{n}- h\;y'(x_{n} ) +\frac{h^2}{2!}\;y''(x_{n}) -\frac{h^3}{3!}\;y'''(x_{n}) +\frac{h^4}{4!} y^{iv}(x_n)+...

Add the two equations together and you get:

y_{n+1} +y_{n-1} =2 y_{n} +h^2 y''(x_{n}) +\frac{h^{4}}{12} y^{iv}(x_n) +O(h^6)

Rearrange:

y_{n+1} -2y_n +y_{n-1} =\mathbf {h^2 y''(x_n)} +\mathbf{\frac{h^4}{12} y^{iv}(x_n)}+O(h^6)

We will largely ignore O(h^6) and substitute for the bold h^2 y''(x_n)
and \frac{h^4}{12} y^{iv}(x_n) above.
For \mathbf{h^2 y''(x_n)} we use the original differential equation:

y''(x_n)=-g_ny_n+s_n

For the bold
\mathbf{\frac{h^4}{12} y^{iv}(x_n)}, we figure that:

\frac{h^4}{12}\frac{d^2}{dx^2} (y''(x_n)) =\frac{h^4}{12}\frac{d^2}{dx^2} (-g_n y_n +s_n)

Now we just learned (rearranging the prior y_{n+1}-2y_n+y_{n-1} equation again), that for generic y(x),

y''(x_n) =\frac{y_{n+1} -2y_n +y_{n-1}-O(h^4)}{h^2}

Use that idea twice to express the pieces of \frac{d^2}{dx^2}(-g_ny_n+s_n) and you get:

\frac{h^4}{12}\frac{d^2}{dx^2} (-g_n y_n +s_n) =\frac{h^4}{12} \left[ \frac{-(g_{n+1} y_{n+1} -2g_n y_n +g_{n-1} y_{n-1})}{h^2} +\frac{s_{n+1} -2s_n +s_{n-1}} {h^2} \right]

…which is then equal to:

\frac{h^2}{12} [-(g_{n+1} y_{n+1} -2g_{n} y_{n} +g_{n-1} y_{n-1}) +s_{n+1} -2s_{n} +s_{n-1}]

Now sub that back into the equation with the boldface:

y_{n+1} -2y_n +y_{n-1} =h^2(-g_n y_n +s_n) +\frac{h^2}{12} [-g_{n+1} y_{n+1} +2g_{n} y_{n} -g_{n-1} y_{n-1} +s_{n+1} -2s_{n} +s_{n-1}]

Divide both sides by h^2:

\frac{y_{n+1} -2y_n +y_{n-1}}{h^2} =(-g_n y_n +s_n) +\frac{1}{12} [-g_{n+1} y_{n+1} +2g_{n} y_{n} -g_{n-1} y_{n-1} +s_{n+1} -2s_{n} +s_{n-1}]

…which is equal to (combining terms):

\frac{-1}{12}(g_{n+1} y_{n+1} +10g_n y_n +g_{n-1} y_{n-1}) +\frac{1}{12} (s_{n+1} +10s_n +s_{n-1})

Now, applying the regular Numerov algorithm might consist of picking two points y_{n-1}
and y_n, using the above formula to generate
y_{n+1} and iterating over all points
n in the x grid. But here I will substitute into the one-dimensional time-independent Schrödinger equation and introduce the matrix method.
Rearrange the Schrödinger Equation:
The one-dimensional Schrödinger equation gives

\frac{-\hbar^2}{2m}\frac{d^2\psi}{dx^2} +V(x)\psi(x)=E\psi(x)

Rearranging,

\frac{d^2\psi}{dx^2}= \frac{2m}{\hbar^2}V(x)\psi(x) -\frac{2mE}{\hbar^2}\psi(x)

Looking familiar yet? Let y''(x)=\frac{d^2\psi}{dx^2},
and -g(x)=\frac{2m}{\hbar^2}V(x),
y(x)=\psi(x),
and s(x)=\frac{-2mE}{\hbar^2}\psi(x).
Now again define a grid of x values (how to get the grid spacing right to come later…) and let \psi_n represent the value of \psi(x)
at the x_n’th point.
Now we can use the result above, noting that \hbar is the reduced Planck constant, but h is the grid spacing \Delta x:

\begin{align} \frac{y_{n+1} -2y_n +y_{n+1}}{h^2} = & \frac{-1}{12} \left[ (\frac{-2m}{\hbar^2}V_{n+1}) (\psi_{n+1} ) +10(\frac{-2m}{\hbar^2} V_n) \psi_n +(\frac{-2m}{\hbar^2}V_{n-1}) (\psi_{n-1} )\right] \\\ & +\frac{1}{12} \left[ \frac{-2mE}{\hbar^2}\psi_{n+1} +10(\frac{-2mE}{\hbar^2})\psi_n +\frac{-2mE}{\hbar^2}\psi_{n-1} \right] \end{align}

Multiply everything by \frac{-\hbar^2}{2m} and rearrange a little:

\frac{-\hbar^2}{2m}\left[ \frac{\psi_{n+1} -2\psi_n +\psi_{n-1}}{h^2} \right] +\frac{1}{12} \left[V_{n+1} \psi_{n+1} +10V_n \psi_n +V_{n-1} \psi_{n-1}\right] =\frac{1}{12}E \left[ \psi_{n+1} +10\psi_n +\psi_{n-1} \right]

And now for…

The Matrix Form

Let \psi(x) be represented by a column vector…
I.e.

\begin{pmatrix} \psi_{1} \\\ \psi_{2} \\\ \psi_{3} \\\ \vdots \end{pmatrix}

Here each \psi_i
is the value of \psi(x)
at the grid point x_i (Many apologies if the above shows up as a row vector; I tried).
Now we can express the equation above in matrix format as…

\frac{-\hbar^2}{2m}\mathbf{A} \;\mathbf{\psi} +\mathbf{B}\;\mathbf{V}\;\mathbf{\psi} =E\;\mathbf{B}\;\mathbf{\psi}

Now, \mathbf{A} is equal to…
\frac{1}{h^2}*

\begin{pmatrix} -2 & 1 & 0 & 0 & 0 &\cdots\\\ 1 & -2 & 1 & 0 & 0 &\cdots\\\ 0 & 1 & -2 & 1 & 0 &\cdots\\\ 0 & 0 & 1 & -2 & 1 &\cdots\\\ 0 & 0 & 0 & 1 & -2 &\cdots\\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix}

And B is \frac{1}{12}*

\begin{pmatrix} 10 & 1 & 0 & 0 & 0 &\cdots\\\ 1 & 10 & 1 & 0 & 0 &\cdots\\\ 0 & 1 & 10 & 1 & 0 &\cdots\\\ 0 & 0 & 1 & 10 & 1 &\cdots\\\ 0 & 0 & 0 & 1 & 10 &\cdots\\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix}

Finally, V is equal to the diagonal matrix:

\begin{pmatrix} V_{1} & 0 & 0 & 0 & 0 &\cdots\\\ 0 & V_{2} & 0 & 0 & 0 &\cdots\\\ 0 & 0 & V_{3} & 0 & 0 &\cdots\\\ 0 & 0 & 0 & V_{4} & 0 &\cdots\\\ 0 & 0 & 0 & 0 & V_{5} &\cdots\\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix}

(Does this remind anyone of Discrete Variable Representation? Mathematica code for that to come soon…)

Left-multiplying by \mathbf{B^{-1}} gives…

\frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A\;\psi} +\mathbf{V\;\psi} =E\;\psi

Rearranging,

\left( \frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A} +\mathbf{V} \right) \psi=E\;\psi

And voila! we have an equation in the \mathbf{H} \psi = E \psi form.
So \frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A} is your kinetic energy operator. If you’re not sure the matrix form is equivalent, pick three elements of the \psi vector to be
\psi_{n-1},
\psi_n,
and \psi_{n+1},
and propagate through the matrix algebra element-wise.

I will now define these matrices and operations in Mathematica…

Mathematica Code

Grid spacing should be fine enough to have five points in the narrowest half-wavelength. For the problem I was working on, I cared about wavefunctions near an energy of 8.8475510^-16 ergs, so p0 = \sqrt{2 m E} = 1.4327*10^-19 →
the narrowest lobe has width \frac{1}{2}*\frac{h}{p0}=2.31244*10^-8 (this time h = Planck’s constant in ergsec).
The grid spacing should be \frac{1}{5} of that width, or 4.62488*10^-9… I believe this problem was using inches for some unfathomable reason.

To use this function, first define your potential function as = , i.e.:

v1input[x_] = a Abs[x]

(Here a was a constant set equal to 5.368244090805358*10^-9, the slope of the V-shaped potential).
Then, decide your minimum and maximum values of x (advice for the V-shaped potential: go a couple wavelengths beyond the wall at the energy of interest…).

Then employ the following function, where expr_ is the name of your potential function…

MatrixNumerov[expr_, xmin_, xmax_, dx_, hbar_, m_] := 
 Module[{xvals, Vfunction, Vmatrix, A, B, KE, n1, H, eivals, eivecs},
  xvals = Table[i, {i, xmin, xmax, dx}];
  Vfunction = expr[xvals];
  Vmatrix = DiagonalMatrix[Vfunction];
  n1 = Length[xvals];
  B = (DiagonalMatrix[1 + 0 Range[n1 - Abs[-1]], -1] + 
      10 DiagonalMatrix[1 + 0 Range[n1 - Abs[0]], 0] + 
      DiagonalMatrix[1 + 0 Range[n1 - Abs[1]], 1])/12.; 
  A = (DiagonalMatrix[1 + 0 Range[n1 - Abs[-1]], -1] - 
      2 DiagonalMatrix[1 + 0 Range[n1 - Abs[0]], 0] + 
      DiagonalMatrix[1 + 0 Range[n1 - Abs[1]], 1])/((dx)^2); 
  KE = (-hbar^2/(2 m))*Inverse[B].A;
  H = KE + Vmatrix;
  {eivals, eivecs} = Eigensystem[H];
  {eivals, eivecs}]

Here’s an example…
With x running from -2.573102897338254`E-7 to 2.573102897338254E-7 in increments of 4.624881629283604’E-9,

the code:

MatrixNumerov[v1input, -2.573102897338254`*^-7, 
2.573102897338254`*^-7, 4.624881629283604`*^-9, 1.0545718`*^-27, 
1.1600000000000001`*^-23]

…returns a list of eigenvalues and a list of eigenvectors, just like a regular Eigensystem[…] command,
with the smallest ten eigenvalues equal to $8.8466910^{-16}, 8.2101310^{-16}, 7.5579610^{-16}, 6.8641110^{-16},
6.1480610^{-16}, 5.3684310^{-16}, 4.552710^{-16}, 3.6180810^{-16},
2.6039610^{-16}, 1.1364110^{-16}$
.
Those ten eigen-energies solved by the Airy method were:
$8.8475510^{-16}, 8.2105510^{-16}, 7.558510^{-16}, 6.864210^{-16},
6.1483610^{-16}, 5.3682410^{-16}, 4.5528310^{-16}, 3.6175810^{-16},
2.60410^{-16}, 1.1346510^{-16}$
So you can see, the approximation is darn close.
This has been adapted from my quantum mechanics midterm (MIT course 5.73 taught by Bob Field Fall 2018) - apologies for any errors/omissions, and I’m afraid I’ll have to leave a discussion of how to get the energies by the Airy functions for V-shaped potential with slope = a (above) for another time.

1 Like
2

Hi its Nice. I have a question, what changes can be done in your’s code to run it for three and four dimensional Potential?

1 Like
3

I’m not sure but when I figure it out I’ll write another article.