Matlabが異なる行列式関数を使用しているため（または2つの異なる方法に含まれる浮動小数点の精度に関係する何かのような他の理由）、小さすぎる行列式の値を探しています。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))
私が実際に行ったことはすべて、オメガンのすべての値の行列式の値をログに記録し、それらの行列式の値のログをオメガンの関数としてプロットしたことです。プロットは次のとおりです。
グラフに3つの大きな落ち込みがあります。2つは16.3818と32.7636の結果と一致しますが、欠落している追加の1つもあります（おそらく、行列式が1e-10未満であるという条件が低すぎて、Fortranコードでそれを取得できなかったためです）。したがって、Matlabは、これらが探していたオメガンの値であることも示していますが、行列式が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)
とオメガのために解決します。これは私がそれについて行った方法です：
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コードで見つけた2つだけでなく、4つすべてのソリューションが得られます（ただし、そのうちの1つであるゼロはomeganのテスト範囲外でした）。これまで行ってきたように、Matlabで行列式をチェックしてこれを解決したい場合は、行列式の絶対値がより小さいことをチェックしている値で遊ぶ必要があります。私はそれを1e-4の値で動作させました（16.382、28.374、および32.764の3つのソリューションを提供しました）。
このような長い解決策については申し訳ありませんが、うまくいけばそれが役立つでしょう。
アップデート：
上記のコードの最初のブロックで、
vals(i) = abs(det(a));
と
[L,U] = lu(a);
s = det(L);
vals(i) = abs(s*prod(diag(U)));
これは、Matlabのドキュメントによるとdetが使用していると思われるアルゴリズムです。これで、条件として1E-10を使用できるようになり、動作します。では、Matlabは、ドキュメントに記載されているとおりに行列式を計算していないのではないでしょうか。これはちょっと気がかりです。