Message Boards Message Boards

GROUPS:

Matrix Numerov Function which can be Generalized to Many Potentials

Posted 2 months ago
461 Views
|
2 Replies
|
2 Total Likes
|

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_n y_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^2) =\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_n y_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.84755*10^-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.6248810^-9... I believe this problem was using inches for some unfathomable reason.

To use this function, first define your potential function as <some name> = <some function>, 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.84669*10^{-16}, 8.21013*10^{-16}, 7.55796*10^{-16}, 6.86411*10^{-16}, 6.14806*10^{-16}, 5.36843*10^{-16}, 4.5527*10^{-16}, 3.61808*10^{-16}, 2.60396*10^{-16}, 1.13641*10^{-16}$ . Those ten eigen-energies solved by the Airy method were: $8.84755*10^{-16}, 8.21055*10^{-16}, 7.5585*10^{-16}, 6.8642*10^{-16}, 6.14836*10^{-16}, 5.36824*10^{-16}, 4.55283*10^{-16}, 3.61758*10^{-16}, 2.604*10^{-16}, 1.13465*10^{-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.

2 Replies
Posted 2 months ago

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?

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

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract