7.4: Numerical Methods - Eigenvalue Problem

For illustrative purposes, we develop our numerical methods for what is perhaps the simplest eigenvalue ode. With \(y=y(x)\) and \(0 \leq x \leq 1\), this simple ode is given by

\[y^{\prime \prime}+\lambda^{2} y=0 \label{7.9} \]

To solve Equation \ref{7.9} numerically, we will develop both a finite difference method and a shooting method. Furthermore, we will show how to solve \((7.9)\) with homogeneous boundary conditions on either the function \(y\) or its derivative \(y^{\prime}\).

Finite difference method

We first consider solving \((7.9)\) with the homogeneous boundary conditions \(y(0)=\) \(y(1)=0 .\) In this case, we have already shown that the eigenvalues of Equation \ref{7.9} are given by \(\lambda=\pi, 2 \pi, 3 \pi, \ldots\).

With \(n\) interior points, we have \(x_{i}=i h\) for \(i=0, \ldots, n+1\), and \(h=1 /(n+1)\). Using the centered-finite-difference approximation for the second derivative, Equation \ref{7.9} becomes

\[y_{i-1}-2 y_{i}+y_{i+1}=-h^{2} \lambda^{2} y_{i} \nonumber \]

Applying the boundary conditions \(y_{0}=y_{n+1}=0\), the first equation with \(i=1\), and the last equation with \(i=n\), are given by

\[\begin{aligned} -2 y_{1}+y_{2} &=-h^{2} \lambda^{2} y_{1} \\ y_{n-1}-2 y_{n} &=-h^{2} \lambda^{2} y_{n} \end{aligned} \nonumber \]

The remaining \(n-2\) equations are given by \((7.10)\) for \(i=2, \ldots, n-1\) It is of interest to see how the solution develops with increasing \(n\). The smallest possible value is \(n=1\), corresponding to a single interior point, and since \(h=1 / 2\) we have

\[-2 y_{1}=-\frac{1}{4} \lambda^{2} y_{1} \nonumber \]

so that \(\lambda^{2}=8\), or \(\lambda=2 \sqrt{2}=2.8284\). This is to be compared to the first eigenvalue which is \(\lambda=\pi\). When \(n=2\), we have \(h=1 / 3\), and the resulting two equations written in matrix form are given by

\[\left(\begin{array}{rr} -2 & 1 \\ 1 & -2 \end{array}\right)\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right)=-\frac{1}{9} \lambda^{2}\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right) \nonumber \]

This is a matrix eigenvalue problem with the eigenvalue given by \(\mu=-\lambda^{2} / 9\). The solution for \(\mu\) is arrived at by solving

\[\operatorname{det}\left(\begin{array}{cc} -2-\mu & 1 \\ 1 & -2-\mu \end{array}\right)=0 \nonumber \]

with resulting quadratic equation

\[(2+\mu)^{2}-1=0 \nonumber \]

The solutions are \(\mu=-1,-3\), and since \(\lambda=3 \sqrt{-\mu}\), we have \(\lambda=3,3 \sqrt{3}=5.1962\). These two eigenvalues serve as rough approximations to the first two eigenvalues \(\pi\) and \(2 \pi\).

With A an \(n\)-by-n matrix, the MATLAB variable mu=eig(A) is a vector containing the \(n\) eigenvalues of the matrix A. The built-in function eig.m can therefore be used to find the eigenvalues. With \(n\) grid points, the smaller eigenvalues will converge more rapidly than the larger ones.

We can also consider boundary conditions on the derivative, or mixed boundary conditions. For example, consider the mixed boundary conditions given by \(y(0)=0\) and \(y^{\prime}(1)=0\). The eigenvalues of Equation \ref{7.9} can then be determined analytically to be \(\lambda_{i}=(i-1 / 2) \pi\), with \(i\) a natural number.

The difficulty we now face is how to implement a boundary condition on the derivative. Our computation of \(y^{\prime \prime}\) uses a second-order method, and we would like the computation of the first derivative to also be second order. The condition \(y^{\prime}(1)=0\) occurs on the right-most boundary, and we can make use of the second-order backward-difference approximation to the derivative that we have previously derived. This finite-difference approximation for \(y^{\prime}(1)\) can be written as

\[y_{n+1}^{\prime}=\frac{3 y_{n+1}-4 y_{n}+y_{n-1}}{2 h} . \nonumber \]

Now, the \(n\)th finite-difference equation was given by

\[y_{n-1}-2 y_{n}+y_{n+1}=-h^{2} y_{n}, \nonumber \]

and we now replace the value \(y_{n+1}\) using \((7.11)\); that is,

\[y_{n+1}=\frac{1}{3}\left(2 h y_{n+1}^{\prime}+4 y_{n}-y_{n-1}\right) . \nonumber \]

Implementing the boundary condition \(y_{n+1}^{\prime}=0\), we have

\[y_{n+1}=\frac{4}{3} y_{n}-\frac{1}{3} y_{n-1} \nonumber \]

Therefore, the \(n\)th finite-difference equation becomes

\[\frac{2}{3} y_{n-1}-\frac{2}{3} y_{n}=-h^{2} \lambda^{2} y_{n} \nonumber \]

For example, when \(n=2\), the finite difference equations become

\[\left(\begin{array}{rr} -2 & 1 \\ \frac{2}{3} & -\frac{2}{3} \end{array}\right)\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right)=-\frac{1}{9} \lambda^{2}\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right) \nonumber \]

The eigenvalues of the matrix are now the solution of

\[(\mu+2)\left(\mu+\frac{2}{3}\right)-\frac{2}{3}=0 \nonumber \]

Or

\[3 \mu^{2}+8 \mu+2=0 \nonumber \]

Therefore, \(\mu=(-4 \pm \sqrt{10}) / 3\), and we find \(\lambda=1.5853,4.6354\), which are approximations to \(\pi / 2\) and \(3 \pi / 2\).

Shooting method

We apply the shooting method to solve Equation \ref{7.9} with boundary conditions \(y(0)=\) \(y(1)=0\). The initial value problem to solve is

\[\begin{aligned} &y^{\prime}=z \\ &z^{\prime}=-\lambda^{2} y \end{aligned} \nonumber \]

with known boundary condition \(y(0)=0\) and an unknown boundary condition on \(y^{\prime}(0)\). In fact, any nonzero boundary condition on \(y^{\prime}(0)\) can be chosen: the differential equation is linear and the boundary conditions are homogeneous, so that if \(y(x)\) is an eigenfunction then so is \(A y(x)\). What we need to find here is the value of \(\lambda\) such that \(y(1)=0 .\) In other words, choosing \(y^{\prime}(0)=1\), we solve

\[F(\lambda)=0 \nonumber \]

where \(F(\lambda)=y(1)\), obtained by solving the initial value problem. Again, an iteration for the roots of \(F(\lambda)\) can be done using the Secant Method. For the eigenvalue problem, there are an infinite number of roots, and the choice of the two initial guesses for \(\lambda\) will then determine to which root the iteration will converge.

For this simple problem, it is possible to write explicitly the equation \(F(\lambda)=0\). The general solution to Equation \ref{7.9} is given by

\[y(x)=A \cos (\lambda x)+B \sin (\lambda x) . \nonumber \]

The initial condition \(y(0)=0\) yields \(A=0 .\) The initial condition \(y^{\prime}(0)=1\) yields

\[B=1 / \lambda \nonumber \]

Therefore, the solution to the initial value problem is

\[y(x)=\frac{\sin (\lambda x)}{\lambda} \nonumber \]

The function \(F(\lambda)=y(1)\) is therefore given by

\[F(\lambda)=\frac{\sin \lambda}{\lambda} \nonumber \]

and the roots occur when \(\lambda=\pi, 2 \pi, \ldots .\)

If the boundary conditions were \(y(0)=0\) and \(y^{\prime}(1)=0\), for example, then we would simply redefine \(F(\lambda)=y^{\prime}(1)\). We would then have

\[F(\lambda)=\frac{\cos \lambda}{\lambda} \nonumber \]

and the roots occur when \(\lambda=\pi / 2,3 \pi / 2, \ldots\)