MATLAB - Eig特徵向量算法

Question

我在仿真中遇到了幾個星期的數值問題，最後我把它縮小到了MATLAB中Eig函數的問題。沒有太多細節，這裏有一個小腳本，它設置了兩個矩陣K1和K，如果你打印它們是完全相同的。但由於某種原因，Eig函數不會爲K返回正確的特徵向量，我不知道爲什麼。在這個矩陣中，我使用了長度爲1的均勻網格上的相鄰點之間的網格距離爲dx = x（i + 1）-x（i），但這與在K1中設置的網格的唯一差別是簡單的使用dx = 1/N。任何人都可以看到地球上發生了什麼？此外，問題似乎只發生在足夠大的N，即小網格距離。我不知道爲什麼，但我對此非常沮喪，因爲它對我的碩士論文至關重要。

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

由於算術運算，形成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'); 

這個額外的代碼生成如下圖所示：