Eigenvectors of an SPD matrix being saved as complex doubles

• I'm creating a reduced-order model for the heat equation in 1-dimension based on the finite element method. Part of the algorithm is producing an N x M matrix, Y where N is the number of spatial nodes and M is the number of time steps being considered.
• Once Y is produced, I compute the eigenvalues and eigenvectors of the matrix Y^t * M * Y, where M is the mass matrix used in the finite element computation.
• The aforementioned matrix, Y^t * M * Y, is ostensibly SPD and should have all real, positive eigenvalues and eigenvectors
• However, upon computing the eig. vals. and eig. vecs. Matlab is saving the values as complex doubles.

I used several functions in this script: phi.m, rhs.m, and simpsons_rule.m, find these files attached

clc
clear
close all
% Initialize variables
tic

Nx = 100;
a = 0;
b = 1;
x = linspace(a,b,Nx);
t0 = 0;
t1 = 0.6;
Nt = 100;
t = linspace(t0,t1,Nt);
dt = (t(2) - t(1));
h = (x(2) - x(1));
Fh = zeros(Nx,1);    % Load vector
Mh = zeros(Nx,Nx);   % Mass matrix
Sh = zeros(Nx,Nx);   % Stiffness matrix
vh = zeros(Nx,1);    % Numerical solution
Y  = zeros(Nx,Nt);   % Snapshot matrix 
%C  = zeros(Nx,R);   
% vr = zeros(Nx,1);
for ii = 1:Nx        % Initialize numerical solution w/ IC
    vh(ii) = ( x(ii)^2 - x(ii) )*sin(x(ii));
end
Y(:,1) = vh;         % Initializing first column. of Y
%% Create the mass matrix
% Numerically integrating phi function using Simpson's rule
xi0 = x(1); xi1 = x(2); xi2 = x(3);
phi_1 = @(x)phi(x,xi0,xi1,xi2);
diag = simpson_rule(phi_1,phi_1,0,1,300);
xf0 = x(2); xf1 = x(3); xf2 = x(4);
phi_2 = @(x)phi(x,xf0,xf1,xf2);
off_diag = simpson_rule(phi_1,phi_2,0,1,300);
for i = 2:Nx-1
    if( i == Nx-1 )
     Mh(i,i) = diag;
    else
    Mh(i,i) = diag;
    Mh(i+1,i) = off_diag;
    Mh(i,i+1) = off_diag;
    end
end

Mh(1,1) = 1;
Mh(Nx,Nx) = 1;

%% Create stiffness matrix
for jj = 2:Nx-2
    Sh(jj,jj) = 2/h;
    Sh(jj,jj+1) = -1/h;
    Sh(jj+1,jj) = -1/h;
end
Sh(Nx-1,Nx-1) = 2/h;
Sh(1,1)= 1;
Sh(Nx,Nx) = 1;

%% Time marching procedure

for k = 1:Nt-1
 rhs1 = @(x) rhs(x,t(k+1));
% Create the load vector 
% Note, the load vector changes in time since it depends on f.
for j = 2:Nx-1
    x0 = x(j-1);
    x1 = x(j);
    x2 = x(j+1);
    phi1 = @(x) phi(x,x0,x1,x2);
   
    Fh(j) = simpson_rule(phi1,rhs1,0,1,300);
end
    vh = (Mh + dt*Sh)\(Mh*vh + dt*Fh);
    Y(:,k+1) = vh;
end
%% ROM Solution
[vita,lambda] = eig(Y'*Mh*Y);  % Eig vecs (vita) and eig. vals (lambda
                               % of  the matrix Y^t Mh Y

Answer 1

[Christine Tobler]

EIG does not recognize the input matrix as symmetric because it's not exactly symmetric. If you compute

A = Y'*Mh*Y
norm(A - A')

you should see some small round-off error there. Use the following to make the input matrix exactly symmetric:

A = Y'*Mh*Y;
A = (A + A')/2;
[vita, lambda] = eig(A);

This will result in real orthgonal eigenvectors, since EIG will now use a specialized algorithm for symmetric matrices.

[Zachary Candelaria]

Beautiful. Thank you so much for the explanation!