Matlab функция определитель пошел наперекосяк

Question

Ниже приводится выдержка из программы шахты:

function [P] = abc(M,f); if det(M) ~= 1, disp(['Matrix M should have determinant 1'])

Я позволяю вариант для пользователя, чтобы не ввести значение «е».

Когда я запускаю abc ([2 1; 1 1]), программа работает нормально и делает то, что она должна. Но когда я запускаю abc ([6 13; 5 11]), мне говорят: «Матрица M должна иметь определитель 1».

Что происходит на Земле?

EDIT: В окне командной строки, я поступил следующим образом:

M = [6 13; 5 11]; if det(M) ~= 1, disp('Im broken'); end

Matlab затем сам сказал мне, что она сломана.

Благодаря

Answer 1

Вы работаете в стандартных проблем, возникающих из-за ограничений, чисел с плавающей точкой. Результатом функции det, вероятно, является что-то вроде 1.000000001.

Общее правило большого пальца: никогда не проверяйте значения с плавающей запятой для равенства.

Answer 2

Добро пожаловать в удивительно дурацкий мир арифметики с плавающей запятой. MATLAB вычисляет определитель, используя разложение LU, т.е. линейную алгебру. Это происходит потому, что детерминант дико неэффективен для массивов даже небольшого размера, если только это не так.

Последствием этого разложения LU является определитель, вычисляемый как число с плавающей запятой. Это не проблема, ЕСЛИ вы вошли в проблему как тривиально простую, как у вас, - определитель матрицы 2x2, состоящий только из небольших целых чисел. В этом случае сам определитель также будет (разумно) малым целым числом. Таким образом, вы можете решить проблему, просто вычислив определитель матрицы 2x2 самостоятельно, используя формулу учебника.

D = A (1,1) * A (2,2) - A (1,2) * A (2,1);

Это будет точно правильно для небольших целых матриц A, хотя даже это может показать некоторую потерю точности для НЕКОТОРЫХ матриц. Например, рассмотрим простой, 2х2 матрица A:

>> A = [1e8 1;1 1e8]; 

Мы знаем, что определитель этой матрицы 1e16-1.

>> det(A) 
ans = 
        1e+16 

Конечно, MATLAB отображает это как 1e16. Но на самом деле число, сгенерированное функцией det в MATLAB, на самом деле составляет 9999999999999998, поэтому 1e16-2. Плохо, если бы я использовал формулу, которую я дал выше для детерминанта 2x2, она вернула бы результат, который по-прежнему неверен, 10000000000000000. Оба результата были отключены на 1. Вы можете узнать больше об этих проблемах, посмотрев на помощь для eps ,

Моя точка зрения, есть некоторые матрицы 2x2, где вычисление определителя будет просто проблематичным, хотя они являются целыми матрицами.

Как только ваши матрицы становятся нецелыми, тогда вещи действительно становятся истинными числами с плавающей запятой, а не целыми числами. Это означает, что вы просто ДОЛЖНЫ использовать сравнения с допусками на них, а не тест для точного единства. В любом случае это хорошее правило.Всегда используйте допуск, когда вы проводите тест на равенство, по крайней мере, пока вы не научитесь достаточно знать, когда не соблюдать это правило!

Таким образом, вы можете выбрать тест, как это:

if abs(det(A) - 1) < (10*eps(1)) 
    warning('The sky is falling! det has failed me.') 
end 

Обратите внимание, что я использовал EPS (1), так как мы сравниваем вещи 1. Тот факт, что я умножил EPS на 10 позволяет wee бит slop в вычислении определителя.

И, наконец, вы должны знать, что независимо от того, какой тест вы используете для определения детерминанта, часто нужно делать BBBBBBBBBBAAAAAAAAAADDDDDDDD! Да, возможно, ваш учитель сказал вам сделать это, или вы нашли что-то в учебнике. Но детерминант - это просто плохая вещь для численных вычислений. Почти всегда альтернативы детерминанту. Опять же, это называется суждением, зная, что то, что вам говорят использовать, на самом деле не так.

Answer 3

Чтобы дать вам представление: det не вычисляется с использованием старой формулы, которую вы изучали в линейной алгебре, но с использованием более эффективных алгоритмов.

Например, используя Gaussian elimination, вы можете преобразовать M в эквивалентную верхнюю треугольную матрицу, а затем вычислить определитель как произведение главной диагонали (являющегося нижним треугольником всех нулей).

M = [6 13; 5 11] 
G = M - [0 0; M(2,1)/M(1,1) * M(1,:)]; 

Теоретически det(M) равно det(G), что 6 * 1/6 = 1, но, будучи Gс плавающей точкой, а не числовой матрицы реального, G(1,1)*G(2,2)~=1!

На самом деле G(1,1) и G(2,2) не точно 1 и 1/6, но они имеют очень малую относительную погрешность (см eps, что на большинстве машин составляет около 2.22e-16). Их реальное значение будет около 6 * (1 + eps) и 1/6 * (1 + eps), поэтому их продукт также будет иметь небольшую ошибку.

Я не уверен, использует ли Matlab исключение Гаусса или аналогичное LU decomposition.