MATLAB - Eig eigenvector algorithm

Question

I have had a numerical problem in my simulations for weeks and I have finally narrowed it down to being a problem with the Eig function in MATLAB. Without going into too many details here is a small script that sets up two matrices K1 and K, which if you print them are completely identical. But for some reason the Eig function does not return the correct eigenvector for K and I have no idea why. In this matrix I have used the grid distance between adjacent points as dx = x(i+1)-x(i) on a uniform grid of length 1 but that is the only difference to the one set up in K1, where I simply used dx=1/N. Can anyone see what on earth is going on? Also the problem seems to happen only for sufficiently large N, i.e. small grid distances.. I have no idea why, but I am so frustrated for this since it is essential for my master's thesis.

    clear all
    clc
    N=1000;
    dx=1/N; %grid distance is 1/N
    x=dx*(1:N); %make grid
    A=zeros(N,N); 
    for i=1:N-1
    A(i,i+1)=1;
    end
    K1=-1/dx^2*(A+A'-2*eye(N));

    for i=2:N-1
    K(i,i)=-2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i))+1/(x(i)-x(i-1)));
    K(i,i-1)=2/(x(i+1)-x(i-1))*(1/(x(i)-x(i-1)));
    K(i,i+1)=2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i)));
    end
    K(N,N-1)=K(N-1,N-2);
    K(1,1)=K(2,2);
    K(1,2)=K(2,3);
    K(2,1)=K(2,3);
    K(N,N)=K(N-1,N-1);
    K=-K;
    [h,y]=eig(K); %eigenvectors (h) and eigenvalues (y) of first K
    [z,v]=eig(K1);   %eigenvectors (z) and eigenvalues of (v) of K1
    plot(x,z(:,1)) %plot first eigenvector of K1
    plot(x,h(:,1))

Answer 1

Since the arithmetic operations to form K and K1 are different, their entries are not necessarily equivalent in floating point arithmetic:

>> norm(K-K1,2)
ans =
   3.8687e-07

Therefore, the eigenvalues will not match exactly nor are they guaranteed to come in the same order, which means the eigenvectors are in different columns:

>> t = [diag(y),diag(v)];
>> t(1:5,:)
ans =
   1.0e+06 *

    3.9190    0.0000
    3.9207    0.0000
    3.9225    0.0001
    3.9242    0.0002
    3.9259    0.0002

However, all of the eigenvalues, when given a similar ordering, have approximately the same value:

>> norm([sort(diag(y))-sort(diag(v))],2)
ans =
   3.4607e-07

>> norm([sort(diag(y))-sort(diag(v))],2)/norm(diag(y),2)
ans =
   4.4685e-15

---

The eigenvalues and eigenvectors of both matrices are, using the metric of relative machine precision, the same. However, the arithmetic operations used to create K1 create a (near-)perfect symmetric matrix which results in smooth eigenvalues from the algorithm used by eig, possibly some Hessenberg transformation since you're requesting eigenvectors as well. K, however, is not symmetric to working precision due to the method of calculation and, if I am reading your intent correctly, will not be for a non-uniform grid.

If you need the eigenvalues to be sorted, you must sort them yourself because general matrices will produce randomly order eigenvalues:

%    Pull diagonals
v = diag(v);
y = diag(y);

%   Sort the non-symmetric eigenvalues in ascending order
[y,orderedIndex] = sort(y);
h = h(:,orderedIndex);

%   Perform a sign correction such that all eigenvectors have the same sense.
%   This is not needed computationally but doing it for easy comparison.
h = abs(h).*sign(z);

%   Plot
subplot(2,1,1);
    plot(1:N,v,1:N,y,'--');
    legend('K1','K','Location','SouthWest');
subplot(2,1,2);
    g = plot(x(:),z(:,1),x(:),h(:,1),'--');
    legend('K1','K','Location','SouthWest');

This additional code generates the following figure: