% Eigenvalues of Differential Equations
% We can use finite difference techniques to solve for the eigenvalues of an ordinary differential equation. We use the finite difference approximation to convert the ODE to a series of linear algebraic equations and then solve for the eigenvalues of the resulting matrix.
% To demonstrate this technique, we will solve the following problem:
%  \begin{array}{l}
% \frac{d^2 y}{dx^2 }+\lambda^2 y=0\\
% y\left(0\right)=y\left(1\right)=0
% \end{array}
% Our goal is to determine the value of \lambda for which there is a solution to this system. For comparison, the analytical solution is
%   \begin{array}{l}
% \lambda_n =n\pi \\
% y_n =A_n \;\sin \left(n\pi x\right)
% \end{array}
% To carry out the numerical solution, we use a finite difference approximation as follows:
% \frac{y_{i-1} -2y_i +y_{i+1} }{h^2 }+\lambda^2 y_i =0
%   which can be rewritten as
%  y_{i-1} -\left(2-\lambda^2 h^2 \right)y_i +y_{i+1} =0
% If we divide the range 0<x<1 into 8 divisions, this yields 9 mesh points (including the end points). The boundary conditions prescribe two of these values, so we are left with a system of 7 equations and 7 unknowns (not including \lambda). The resulting system can be written as
%   \left\lbrack \begin{array}{ccccccc}
% -2 & 1 & 0 & 0 & 0 & 0 & 0\\
% 1 & -2 & 1 & 0 & 0 & 0 & 0\\
% 0 & 1 & -2 & 1 & 0 & 0 & 0\\
% 0 & 0 & 1 & -2 & 1 & 0 & 0\\
% 0 & 0 & 0 & 1 & -2 & 1 & 0\\
% 0 & 0 & 0 & 0 & 1 & -2 & 1\\
% 0 & 0 & 0 & 0 & 0 & 1 & -2
% \end{array}\right\rbrack \left\lbrace y\right\rbrace +\lambda^2 h^2 \left\lbrace y\right\rbrace =0
% Where \left\lbrace y\right\rbrace is a vector of y-values at the mesh points and, given that we used 8 divisions, h=1/8.
% If we define the variable e_i as the list of eigenvalues of the 7x7 matrix in the previous equation, then we can write (based on the definition of the eigenvalues of a matrix)
% \lambda_i^2 h^2 =-e_i
% Our approach, then, is to use the built-in eig command to find the eigenvalues of the matrix in the equation above and then find the values of \lambda corresponding to each value of e_i. In this case, we’ll get 7 values for \lambda. If we use more divisions on the interval 0<x<1, then we will get more values for \lambda.
% Sample code for this problem is as follows:

h=1/8;
mymat=[-2 1 0 0 0 0 0;
       1 -2 1 0 0 0 0;
       0 1 -2 1 0 0 0;
       0 0 1 -2 1 0 0;
       0 0 0 1 -2 1 0;
       0 0 0 0 1 -2 1;
       0 0 0 0 0 1 -2];
%Below we calculate the eigenvalues and sort them for convenience
D=sort(eig(mymat),'descend');
lam=sqrt(-D)/h;
check=lam/pi;
myint=(1:7)';
plot(myint,check,myint,myint, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('normalized eigenvalues - \lambda /pi')
myerr=abs(check-myint)./(myint);
semilogy(myint,myerr, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('relative error')

% You can see that the lowest eigenvalues are the most accurate and the error gets quite large for the larger ones.
% Note that we can write this to allow for a variable number of divisions if we are careful about how we define the matrix. Increasing the number of divisions gives us more eigenvalues and better accuracy.

n=7;
h=1/(n+1);
voffdiag=ones(n-1,1);
mymat=-2*eye(n)+diag(voffdiag,1)+diag(voffdiag,-1);
D=sort(eig(mymat),'descend');
lam=sqrt(-D)/h;
check=lam/pi;
myint=(1:n)';
plot(myint,check,myint,myint, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('normalized eigenvalues - \lambda /pi')
myerr=abs(check-myint)./(myint);
semilogy(myint,myerr, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('relative error')

%  We see reduced error for the smaller eigenvalues as we increase the number of divisions.

n=100;
h=1/(n+1);
voffdiag=ones(n-1,1);
mymat=-2*eye(n)+diag(voffdiag,1)+diag(voffdiag,-1);
D=sort(eig(mymat),'descend');
lam=sqrt(-D)/h;
check=lam/pi;
myint=(1:n)';
plot(myint,check,myint,myint, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('normalized eigenvalues - \lambda /pi')
myerr=abs(check-myint)./(myint);
semilogy(myint,myerr, 'LineWidth',2)
xlabel('eigenvalue number')
ylabel('relative error')

% Note that the code is written such that the number of points can easily be varied.
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu