我在仿真中遇到了幾個星期的數值問題,最後我把它縮小到了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'之間的區別。當然,一般來說,它們會有不同的特徵值 –