Вы ищете слишком малые значения детерминанта, потому что 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))
Все, что я на самом деле сделал, - это записал значения определителя для всех значений омегана и составил график этих значений определителя как функции омегана. Вот сюжет:
Вы заметили три основных провала на графике. Два совпадают с вашими результатами 16.3818 и 32.7636, но есть и дополнительный, который вы пропустили (возможно, потому что ваше условие определителя было меньше 1e-10, было слишком низким даже для вашего кода на Фортране, чтобы его забрать). Поэтому 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
Это дает вам все четыре решения, а не только два, которые вы нашли с вашим кодом на Фортране (хотя одно из них, ноль, было за пределами диапазона тестирования для 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)));
- это алгоритм, который det предположительно использует в соответствии с документами Matlab. Теперь я могу использовать 1E-10 в качестве условия, и оно работает. Так что, может быть, Matlab не рассчитывает определитель точно так, как говорят документы? Это немного тревожно.