一千萬個為什麽

搜索

matlab精度決定因素問題

我有以下程序

format compact; format short g; clear; clc;  
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;  
for i=1:500000  
omegan=1.+0.0001*i;

a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0;
a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0;
a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1;
a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2;

if(abs(det(a))<1E-10) sprintf('omegan= %8.3f det= %8.3f',omegan,det(a))
end          
end

上述系統的分析解決方案,以及同樣的用fortran編寫的程序給出了omegan的值等於16.3818和32.7636(fortran值;分析略有不同,但它們在某處)。

那麽,現在我想知道......我哪裏出錯了?為什麽matlab沒有給出預期的結果?

(這可能非常簡單,但令我頭疼)

最佳答案

您正在尋找太小的行列式值,因為Matlab使用的是不同的行列式函數(或者某些其他原因,例如與兩種不同方法中涉及的浮點精度有關)。我將向您展示Matlab實際上為您提供了正確的值,並且更好地解決了這個問題。

首先,讓我們采取您的代碼並稍微改變它。

format compact; format short g; clear; clc;
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
vals = zeros(1,500000);
for i=1:500000
    omegan=1.+0.0001*i;

    a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0;
    a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0;
    a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1;
    a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2;
    vals(i) = abs(det(a));
    if(vals(i)<1E-10)
        sprintf('omegan= %8.3f det= %8.3f',omegan,det(a))
    end
end
plot(1.+0.0001*(1:500000),log(vals))

我所做的一切都是記錄了所有omegan值的行列式值,並將這些決定因素值的對數繪制為omegan的函數。這是情節:

alt text

您註意到圖中有三個主要下降。兩個與16.3818和32.7636的結果一致,但是還有一個你缺少的(可能是因為你的行列式小於1e-10的條件太低了,即使你的Fortran代碼拿起它也是如此)。因此,Matlab也告訴你那些是你正在尋找的omegan的值,但是因為決定因素是在Matlab中以不同的方式確定的,所以這些值並不相同 - 在處理條件差的矩陣時也就不足為奇了。此外,正如其他人所說,它可能與使用單精度浮標的Fortran有關。我不打算研究為什麽他們不是因為我不想浪費我的時間。相反,讓我們看看你想要做什麽,並嘗試不同的方法。

正如我相信你所知,你正試圖找到矩陣的特征值

a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0 0 2 -2]];

,把它們等於

-omegan^2*(Jm/(G*J)*d^2)

並解決omegan。這就是我的方式:

format compact; format short g; clear; clc;
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
C1 = (Jm/(G*J)*d^2);
a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0,0,2,-2]];
myeigs = eig(a);
myeigs(abs(myeigs) < eps) = 0.0;
for i=1:4
    sprintf('omegan= %8.3f', sqrt(-myeigs(i)/C1))
end

這為您提供了所有四種解決方案 - 不僅僅是您使用Fortran代碼找到的兩種解決方案(盡管其中一種,零,超出了omegan的測試範圍)。如果你想通過檢查Matlab中的行列式來解決這個問題,就像你一直試圖做的那樣,那麽你將不得不使用你正在檢查行列式的絕對值小於的值。我得到它的值為1e-4(它給出了3個解決方案:16.382,28.374和32.764)。

很抱歉這麽長的解決方案,但希望它有所幫助。

更新</強>

在我上面的第一個代碼塊中,我替換了

vals(i) = abs(det(a));

[L,U] = lu(a);
s = det(L);
vals(i) = abs(s*prod(diag(U)));

這是根據Matlab文檔使用det的算法。現在,我可以使用1E-10作為條件,它的工作原理。那麽也許Matlab不像文檔所說的那樣完全計算行列式?這有點令人不安。

轉載註明原文: matlab精度決定因素問題