Вы ищете слишком маленькие значения определителя, потому что 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 было слишком низким даже для вашего кода Fortran, чтобы подобрать его). Следовательно, Matlab также сообщает вам, что это значения омегана, которые вы искали, но из-за того, что определитель был определен в Matlab по-другому, значения не были одинаковыми, что неудивительно при работе с плохо обусловленными матрицами. . Кроме того, это, вероятно, связано с тем, что Фортран использует поплавки одинарной точности, как сказал кто-то другой. Я не собираюсь разбираться, почему это не так, потому что не хочу тратить на это время. Вместо этого давайте посмотрим, что вы пытаетесь сделать, и попробуем другой подход.
Вы, как я уверен, в курсе, пытаетесь найти собственные значения матрицы
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
Это дает вам все четыре решения, а не только два, которые вы нашли с помощью своего кода на Фортране (хотя одно из них, нулевое, находилось за пределами вашего диапазона тестирования для омеган). Если вы хотите решить эту проблему, проверив определитель в 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 не вычисляет определитель точно так, как говорят документы? Это немного беспокоит.
person
Justin Peel
schedule
07.05.2010
det(a)может быть почти таким же, пока оно выше, скажем, 1E-3, и что Fortran систематически возвращает значения ближе к нулю для определителя, если он становится маленьким. Это может дать вам представление о том, как вам нужно адаптировать допуски при переносе из одной программы в другую. - person Jonas schedule 07.05.2010