определяющая проблема точности matlab

Новый ответ:

Вы можете исследовать эту проблему с помощью символьных уравнений, которые дадут мне правильные ответы:

>> clear all             %# Clear all existing variables
>> format long           %# Display more digits of precision
>> syms Jm d omegan G J  %# Your symbolic variables
>> a = ((Jm*(d*omegan)^2)/(G*J)-2).*eye(4)+...  %# Create the matrix a
       diag([2 1 1],1)+...
       diag([1 1 2],-1);
>> solns = solve(det(a),'omegan')  %# Solve for where the determinant is 0

solns =

                                0
                                0
            (G*J*Jm)^(1/2)/(Jm*d)
           -(G*J*Jm)^(1/2)/(Jm*d)
       -(2*(G*J*Jm)^(1/2))/(Jm*d)
        (2*(G*J*Jm)^(1/2))/(Jm*d)
  (3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d)
 -(3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d)

>> solns = subs(solns,{G,J,Jm,d},{8e7,77,10540,140/3})  %# Substitute values

solns =

                   0
                   0
  16.381862247021893
 -16.381862247021893
 -32.763724494043785
  32.763724494043785
  28.374217734436371
 -28.374217734436371

Я думаю, что вы либо просто не знали » t выбор значений в вашем цикле, достаточно близких к решениям для omegan , или ваш порог того, насколько близок определитель к нулю, является слишком строгим. Когда я вставляю указанные значения в a вместе с omegan = 16.3819 (что является ближайшим значением к одному решению, которое дает ваш цикл), я получаю следующее:

>> det(subs(a,{omegan,G,J,Jm,d},{16.3819,8e7,77,10540,140/3}))

ans =

    2.765476845475786e-005

еще больше по абсолютной амплитуде, чем 1e-10.

Вы ищете слишком малые значения детерминанта, потому что 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

Это дает вам все четыре решения, а не только два, которые вы нашли с помощью своего кода на Фортране (хотя одно из них, нулевое, находилось за пределами вашего диапазона тестирования для 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 не вычисляет определитель точно так, как говорится в документации? Это немного тревожно.

Я поставил это как ответ, потому что не могу вставить это в комментарий: Вот как Matlab вычисляет определитель . Я предполагаю, что ошибки округления возникают из-за вычисления произведения нескольких диагональных элементов в U.

Алгоритм

Определитель вычисляется из треугольных множителей, полученных с помощью исключения Гаусса

[L,U] = lu(A) s =  det(L)        
%# This is always +1 or -1  
det(A) = s*prod(diag(U))