您正在寻找的行列式值太小，因为 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 的函数。这是情节：
您注意到图表中的三个主要下降。两个与您的 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)
并解决欧米茄。这就是我的做法：
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没有像文档所说的那样计算行列式？这有点令人不安。