Eigenvectors of an SPD matrix being saved as complex doubles
1 view (last 30 days)
Show older comments
Zachary Candelaria
on 21 Mar 2020
Commented: Zachary Candelaria
on 23 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));
% 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
0 Comments
Accepted Answer
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)
See Also
Categories
Find more on Language Fundamentals in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!