Самое быстрое численное решение вещественного кубического полинома?

Question

R вопрос: В поисках самого быстрого способа ЧИСЛИЧНО решите группу произвольных кубик, у которых, как известно, есть действительные коэффициенты и три действительных корня. Сообщается, что функция polyroot в R использует алгоритм 419 Дженкинса-Трауба для сложных полиномов, но для реальных полиномов авторы ссылаются на свои более ранние работы. Каковы более быстрые варианты для реального кубика или, в более общем случае, для реального полинома?

Answer 1

Численное решение для выполнения этого многократно надежным и стабильным образом включает: (1) Формируют сопутствующую матрицу, (2) находят собственные значения сопутствующей матрицы.

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

Для многочлена:

p(t) = c0 + c1 * t + c2 * t^2 + t^3

сопутствующая матрица:

[[0 0 -c0],[1 0 -c1],[0 1 -c2]]

Найти собственные значения такой матрицы; они соответствуют корням исходного многочлена.

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

Обратите внимание, что коэффициент t^3 равен единице, если это не так в ваших полиномах, вам придется разделить все это на коэффициент и затем продолжить.

Удачи.

Редактировать: Numpy и октава также зависят от этой методологии для вычисления корней полиномов. Смотрите, например, эту ссылку .

Answer 2

Самый быстрый из известных (как мне известно) способов найти реальные решения системы произвольных многочленов от n переменных - это многогранная гомотопия. Подробное объяснение, возможно, выходит за рамки ответа StackOverflow, но по сути это алгоритм пути, который использует структуру каждого уравнения с использованием торических геометрий. Google выдаст вам количество документов .

Возможно, этот вопрос лучше подходит для mathoverflow ?

Answer 3

Изложите ответ Ариетты выше:

> a <- c(1,3,-4)
> m <- matrix(c(0,0,-a[1],1,0,-a[2],0,1,-a[3]), byrow=T, nrow=3)
> roots <- eigen(m, symm=F, only.values=T)$values

Является ли это быстрее или медленнее, чем использование кубического решателя в пакете GSL (как предложено выше knguyen), является вопросом его тестирования в вашей системе.

Answer 4

1) Решите для производного полинома P ', чтобы найти три корня. Смотрите там , чтобы знать, как это сделать правильно. Назовите эти корни a и b (с 2) Для среднего корня используйте несколько шагов деления пополам между a и b, и, когда вы окажетесь достаточно близко, закончите методом Ньютона.

3) Для минимального и максимального корня, «охота» решение. Для максимального корня:

• Начните с x0 = b, x1 = b + (b - a) * лямбда, где лямбда - это умеренное число (скажем, 1.6)
• do x_n = b + (x_ {n - 1} - a) * лямбда, пока P (x_n) и P (b) не будут иметь разные знаки
• Выполнить деление пополам + ньютон между x_ {n - 1} и x_n

Answer 5

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

Answer 6

Вы пытались просмотреть пакет GSL http://cran.r -project.org / web / packages / gsl / index.html ?

Answer 7

Доступны общие методы: метод Ньютона, метод деления пополам, секущий, итерация с фиксированной точкой и т. Д. Google любой из них.

Если у вас есть нелинейная система с другой стороны (например, система по N полиномиальным уравнениям в N неизвестных), такой метод, как старший ньютон высокого порядка б.