For a complete explanation of what this program is doing, please see the related paper
in the American Journal of Physics
(preprint at arXiv:1905.13269 [physics.ed-ph]).
Here is a brief summary \(\ldots\)
Our goal is to solve the time-independent Schrödinger equation (TISE),
\begin{equation}
\hat{H} \psi_n(x) = E_n \psi_n(x),
\end{equation}
where the Hamiltonian operator \(\hat{H}\)\( = -\frac{1}{2} \frac{d^2}{dx^2} + V(x)\) (in units where \(m = ℏ = 1\)) includes the potential energy function
shown above (right).
To do so we expand \(\psi_n(x)\) in terms of a set of normalized sine-wave basis functions, also initially shown above.
Denoting these functions as \(\varphi_m(x)\), we write the expansion as
\begin{equation}
\psi_n = \sum_m c_{mn} \varphi_m,
\end{equation}
where the first subscript on the coefficient \(c_{mn}\) indicates which component, while the second subscript
indicates which \(\psi_n\) eigenfunction. Inserting this expansion on both sides of the TISE and then taking
the inner product with an arbitrary basis function \(\varphi_l\), we obtain (after a few algebraic steps) the TISE in matrix form:
\begin{equation}
\begin{bmatrix}
H_{11} & H_{12} & \cdots \\
H_{21} & H_{22} & \cdots \\
\vdots & \vdots & \ddots \\
\end{bmatrix}
\begin{bmatrix}
c_{1n} \\
c_{2n} \\
\vdots
\end{bmatrix}
= E_n
\begin{bmatrix}
c_{1n} \\
c_{2n} \\
\vdots
\end{bmatrix}.
\end{equation}
Here the matrix elements \(H_{mn}\) are defined by integrals,
\begin{equation}
H_{mn} = \langle \varphi_m | \hat{H} | \varphi_n \rangle = \int\!\varphi_m(x)\hat{H}\varphi_n(x)\,dx,
\end{equation}
which the program computes when it loads, and whenever the matrix is reset. The results are displayed
above (left) as colored squares.
To solve the matrix eigenvalue problem we use a method
devised long ago by Jacobi. Note that the normalized column eigenvectors form an
orthogonal matrix \(C\), which we break down into a product of elementary rotation matrices \(R_i\):
\begin{equation}
C = \begin{bmatrix}
c_{11} & c_{12} & \cdots \\
c_{21} & c_{22} & \cdots \\
\vdots & \vdots & \ddots \\
\end{bmatrix} = R_1 R_2 R_3\ldots,
\end{equation}
where the subscripts on the \(R_i\) symbols indicate the sequence of rotations. Each rotation matrix has the form
\begin{equation}
R = \begin{bmatrix}
\,1 \\[-5pt]
& \ddots \\
&& \cos\theta & \cdots & {-}\sin\theta \\[-3pt]
&& \vdots & 1 & \vdots \\
&& \sin\theta & \cdots & \ \cos\theta \\[-3pt]
&&&&& \ddots \\
&&&&&& 1\, \\
\end{bmatrix},
\end{equation}
consisting of an identity matrix except for the four entries that accomplish a rotation by \(\theta\)
in a particular plane. When you click on an off-diagonal matrix cell above, you are selecting the
plane for one of these rotations. Turning the dial changes the rotation angle \(\theta\).
Each rotation modifies the Hamiltonian matrix according to
\begin{equation}
H_\textrm{new} = R^T H_\textrm{old} R,
\end{equation}
so if you choose the rotation angle to (approximately) zero out an off-diagonal matrix element,
the new Hamiltonian will be closer to diagonal than the old one. Eventually, after enough of these
rotations, you will have diagonalized the Hamiltonian so that its diagonal elements are the energy
eigenvalues.
Meanwhile, the program accumulates the product of rotation matrices to form successive approximations
to the matrix \(C\) of eigenvectors, and uses these (approximate) eigenvectors to construct and
plot the (approximate) eigenfunctions \(\psi_n(x)\) as the calculation proceeds.
The Rotate button calculates the optimum rotation angle according to the formula
\begin{equation}
\theta = \frac{1}{2} \tan^{-1} \biggl( \frac{2H_{mn}}{H_{mm} - H_{nn}} \biggr),
\end{equation}
where \(m\) and \(n\) are the chosen row and column numbers. To derive this formula, multiply out
the matrix product \(R^T H R\) in the rotating two-dimensional subspace, and set either off-diagonal element to zero.