Eigenvectors of an SPD matrix being saved as complex doubles

1 view (last 30 days)
  • 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

Accepted Answer

Christine Tobler
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.

More Answers (0)

Categories

Find more on Language Fundamentals in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!