Как правильно найти полиномиальные корни?

Question

Рассмотрим многочлен, такие как:

p = [1 -9 27 -27]; 

очевидно вещественный корень 3:

polyval(p,3) 

0 

При использовании функции roots

q = roots([1 -9 27 -27]); 

с format short:

q = 

    3.0000 + 0.0000i 
    3.0000 + 0.0000i 
    3.0000 - 0.0000i 

и проверить, если корни вещественны:

bsxfun(@eq,ones(size(q)),isreal(q)) 

0 
0 
0 

И еще хуже с format long я получаю:

roots([1 -9 27 -27]) 

ans = 

    3.000019414068325 + 0.000000000000000i 
    2.999990292965843 + 0.000016813349886i 
    2.999990292965843 - 0.000016813349886i 

Как я могу вычислить корни полинома правильно?

Answer 1

Это связано с неточностями с плавающей запятой. Посмотрите на этот пост для деталей: Is floating point math broken?

Одна вещь, которую вы можете сделать, это округлить ответ/с Шифрование до некоторых знаков после запятой, как это:

q = round(roots([1 -9 27 -27]), 4) % rounding off to 4 decimal places 

Answer 2

Вы, возможно, придется работать символически. Для этого вам нужен Symbolic Math Toolbox.

1. Определить полином как символическую функцию. Вы можете (a) использовать poly2sym для генерации символьного многочлена из его коэффициентов. Или (б) еще лучше, определите символическую функцию напрямую, используя строку. Таким образом, вы избегаете потери точности, которая может возникнуть в результате представления коэффициентов как double.

2. Использовать solve, что символически решает алгебраические уравнения.

код с опцией (а):

p = [1 -9 27 -27]; 
ps = poly2sym(p); 
rs = solve(ps); 

кодекса с опцией (б):

ps = sym('x^3-9*x^2+27*x-27'); 
rs = solve(ps); 

В любом случае, результат является символическим:

>> rs 
rs = 
3 
3 
3 

Возможно, вы захотите преобразовать в числовые значения, используя

r = double(rs); 

В вашем примере, это дает

>> format long 
>> r 
r = 
    3 
    3 
    3 

Answer 3

Это очень специфичное для вашего полинома. В общем, вы должны ожидать, что корень кратности m имеет относительную ошибку с плавающей запятой величиной mu^(1/m), где mu=1e-15 - точность машины.В этом случае кратность равна m=3, поэтому ошибка в диапазоне 10^(-5). Это точно масштаб ошибок в ваших результатах.

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

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

---

Математически, у вас есть некоторый многочлен

p(x)=(x-a)^m*q(x) 

с корнем в x=a кратности m. Из-за операции с плавающей запятой, решатель эффективно «видит» полиномиальное

p(x)+e(x) 

где коэффициенты e(x) имеют размер, который величина коэффициентов p раз mu. Близко к корню a, это возмущенные полином можно эффективно заменить на

(x-a)^m*q(a)+e(a) = 0 <==> (x-a)^m = -e(a)/q(a) 

так, что решения образуют м однонаправленным правильный многоугольник или звезду с центром в точке с радиусом a|e(a)/q(a)|^(1/m), который должен находиться в области |a|*mu^(1/m).