# Eigenvectors of an SPD matrix being saved as complex doubles

5 views (last 30 days)
Zachary Candelaria on 21 Mar 2020
• 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));
% 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

Christine Tobler on 23 Mar 2020
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.

#### 1 Comment

Zachary Candelaria on 23 Mar 2020
Beautiful. Thank you so much for the explanation!