Решение квартического уравнения

10

Существует ли открытая C-реализация для решения квартичных уравнений:

ax+bx³+cx²+dx+e=0

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

Мне нужны только реальные решения.

highsciguy
источник
Вам нужны все (реальные) решения одновременно? Как говорит GertVdE ниже, если у вас есть проблемы со стабильностью решения с закрытой формой, на самом деле нет веской причины не использовать какой-либо алгоритм поиска корня.
Годрик Провидец
3
Я нахожу забавным, что это было помечено как нелинейная алгебра, поскольку вы можете просто вычислить собственные значения матрицы-компаньона, которая уже находится в форме Гессенберга, и применение разверток QR будет довольно простым.
Виктор Лю
2
Взгляните на решатели кубических / квартальных чисел, опубликованные в ACM TOMS (алгоритм 954) . Код, который попадает в этот журнал, обычно очень высокого качества. Сам документ находится за платным доступом, но код можно скачать по этой ссылке .
GoHokies
... (позже отредактируйте) код ACM написан на FORTRAN 90, но мое первое впечатление состоит в том, что его можно вызвать из C, не прилагая больших усилий.
GoHokies
1
@ GoHokies Я думаю, вы должны преобразовать свой комментарий в ответ, потому что я думаю, что это хороший ответ на этот вопрос. Тем более что связанная статья позволяет избежать обычной численной нестабильности, и это абсолютно не тривиальная вещь.
Кирилл

Ответы:

20

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

Классическим примером является пример для квадратного уравнения . Вычисляя корни как вы попадете в беду для многочленов, где с тех пор вы получаете отмену в числитель. Вам нужно вычислить .ax2+bx+c=0

x1,2=b±b24ac2a
b4ac
x1=(b+sign(b)b24ac)2a;x2=ca1x1

Хаймам в своем шедевре «Точность и стабильность численных алгоритмов» (2-е изд, SIAM) использует метод прямого поиска, чтобы найти коэффициенты кубического полинома, для которого классическое аналитическое кубическое решение дает очень неточные результаты. В качестве примера он приводит . Для этого многочлена корни хорошо отделены, и, следовательно, проблема не является плохо обусловленной. Однако, если он вычисляет корни, используя аналитический подход, и оценивает многочлен от этих корней, он получает остаток от , используя метод стабильного стандарта (метод матрицы-компаньона) остаток имеет порядок[a,b,c]=[1.732,1,1.2704]O(102)O(1015), Он предлагает небольшую модификацию алгоритма, но даже тогда он может найти набор коэффициентов, ведущих к вычетам что определенно не годится. См. P480-481 вышеупомянутой книги.O(1011)

В вашем случае я бы применил метод Барстоу . Он использует итеративную комбинацию итераций Ньютона на квадратичных формах (а затем решаются корни квадратичных) и дефляцию. Это легко реализовать, и в Интернете даже есть несколько реализаций.

GertVdE
источник
1
Не могли бы вы объяснить, что вы подразумеваете под «Я бы настоятельно рекомендовал не использовать закрытые формы решений, поскольку они, как правило, очень нестабильны в числовом выражении». Это относится только к полиномам 4-й степени или это общее правило?
NoChance
@ EmmadKareem Я обновил свой ответ выше.
GertVdE
3

Смотрите эти:

LHF
источник
2
Используя этот код для многочлена с коэффициентами, приведенными в моем ответе, я нахожу следующее: , который имеет относительную ошибку по сравнению с реальным корнем (рассчитывается с использованием команды корней Octave, которая использует метод сопутствующей матрицы). Он имеет остаток то время как корень метода сопутствующей матрицы имеет остаток . До вас ли это достаточно хорошо (для компьютерной графики это может быть, для некоторых других приложений, это не будет)x1=1.602644912244132e+00O(108)O(107)O(1015)
GertVdE
1

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

Nemocopperfield
источник
Я только что получил корень кубического примера, процитированного с точностью до 2e-16 (немного выше точности моих операций), используя численные рецепты в кубических формулах c (press et al). Так что есть основания надеяться.
Немокопперфильд