MATLAB - Eig特征向量算法

问题描述：

我在仿真中遇到了几个星期的数值问题，最后我把它缩小到了MATLAB中Eig函数的问题。没有太多细节，这里有一个小脚本，它设置了两个矩阵K1和K，如果你打印它们是完全相同的。但由于某种原因，Eig函数不会为K返回正确的特征向量，我不知道为什么。在这个矩阵中，我使用了长度为1的均匀网格上的相邻点之间的网格距离为dx = x（i + 1）-x（i），但这与在K1中设置的网格的唯一差别是简单的使用dx = 1/N。任何人都可以看到地球上发生了什么？此外，问题似乎只发生在足够大的N，即小网格距离。我不知道为什么，但我对此非常沮丧，因为它对我的硕士论文至关重要。MATLAB - Eig特征向量算法

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))

你能否更好地描述'K'和'K1'之间的区别。当然，一般来说，它们会有不同的特征值 –

答

由于算术运算，形成K和K1是不同的，他们的作品都没有在浮点运算一定等同：

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

因此，特征值将不完全匹配，也不是保证按照相同的顺序排列，这意味着特征向量在不同的列中：

>> 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

然而，所有的特征值，给出一个类似的排序时，具有大致相同的值：

>> 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

特征值和两个矩阵是的本征向量，使用相对机器精度的度量中，相同的。然而，用于创建K1的算术运算创建了一个（接近）完美对称矩阵，由于您要求特征向量，因此可以使用由eig使用的算法生成的平滑特征值，可能还有一些Hessenberg变换。然而，由于计算方法的原因，K与工作精度并不对称，如果我正确读取您的意图，将不适用于非均匀网格。

如果您需要的特征值进行排序，你必须自己对它们进行排序，因为一般的矩阵会产生随机排序的特征值：

% 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');

这个额外的代码生成如下图所示：

MATLAB - Eig特征向量算法

相关推荐