How to numerically obtain the eigenvalues of the second-order derivative from the discretization matrix

eigenvalues-eigenvectorslinear algebraMATLABmatricestridiagonal-matrices

I want to find the eigenvalues of the second-order derivative (with pure Neumann boundary conditions, on the interval $[0,1]$) numerically and compare them with the eigenvalues derived analytically.

Analytical eigenvalues for this problem in the continuous case, according to this Wiki page, are the following:

$$\lambda_j = -(j-1)^2 \pi^2,~~~j=1,2,3,…$$

Now let's define the discretized second-order derivative with Neumann boundary conditions in a matrix form. To do this, I basically discretize the equation $\frac{\partial^2 y}{\partial x^2}=0$, writing it as $\mathbf A y_i = 0$, where

$$\mathbf A =
\begin{bmatrix}
-1 & 1 & 0 & 0 & \cdots & 0 \\
1 & -2 & 1 & 0 & \cdots & 0 \\
0 & 1 & -2 & 1 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & \\
0 & \cdots & 0 & 1 & -2 & 1 \\
0 & \cdots & 0 & 0 & -1 & 1 \\
\end{bmatrix}$$

If I try to find the eigenvalues of this matrix numerically, for example, in MATLAB, using the built-in function eig(), I get the following eigenvalues:

Numerical eigenvalues

As you can see, they are not even close to the $\lambda_j$ values for the continuous case. I feel like I am doing something fundamentally wrong here. How can I obtain $\lambda_j$ numerically from the discretization matrix?

Why are eigenvalues of the discretized second-order derivative (written in the matrix form given above) different from the analytical eigenvalues in the continuous case?

The MATLAB code to produce the figure above, just in case:

N = 100;
A = sparse(N, N);

A(1, 1) = -1.0;
A(1, 2) = 1.0;
for i = 2:N-1
    A(i, i-1) = 1.0;
    A(i, i) = -2.0;
    A(i, i+1) = 1.0;
end
A(N, N-1) = -1.0;
A(N, N) = 1.0;

eigvals = eig(full(A));
eigvals = sort(eigvals);

figure('Name', 'Eigenvalues')
plot(eigvals, '.', 'MarkerSize', 10)
grid on

Best Answer

Let me answer my own question. It turns out pretty simple.

First, in the discretization matrix $\mathbf A$, the first and last rows should be multiplied by $h$ (the grid step size), which is equal to $1 / (N - 1)$, where $N$ is the number of grid points. This comes from the fact the first and last rows represent the first derivative, whereas the rest of the rows represent the second derivative, and we multiply the matrix $\mathbf A$ by $h^2$.

Second, since we multiplied the discretization matrix by $h^2$, the eigenvalues will be multiplied by $h^2$ too: $\mathbf A \mathbf y = h^2 \lambda \mathbf y$. This means we need to compare our numerical solution with

$$\lambda_j = -(j-1)^2 \pi^2 h^2,~~~j=1,2,3,...$$

Finally, since the eigenvalues are negative, in the MATLAB code, we need to sort them in descending order.

After fixing these three points, the first (smallest in magnitude) numerical eigenvalues will match the corresponding analytic eigenvalues:

Updated numerical eigenvalues vs analytical solution

Notably, if we increase the number of mesh points $N$, we will get a more accurate solution.

Here is the updated version of the MATLAB code:

N = 100;
h = 1.0 / (N - 1);
A = sparse(N, N);

A(1, 1) = -1.0 * h;
A(1, 2) = 1.0 * h;
for i = 2:N-1
    A(i, i-1) = 1.0;
    A(i, i) = -2.0;
    A(i, i+1) = 1.0;
end
A(N, N-1) = -1.0 * h;
A(N, N) = 1.0 * h;

eigvals = eig(full(A));
eigvals = sort(eigvals, "descend");

figure('Name', 'Eigenvalues')

j = 1:N;
plot(-((j - 1) * pi).^2 * h^2, 'xr', 'MarkerSize', 5)
hold on
plot(eigvals, '.b', 'MarkerSize', 10)
grid on
Related Question