Корни четвертичной функции

Question

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

Я написал функцию, которая, кажется, работает нормально, используя общее решение Ferrari, как показано здесь: http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution.

Вот моя функция:

    private function solveQuartic(A:Number, B:Number, C:Number, D:Number, E:Number):Array{          
        // For paramters: Ax^4 + Bx^3 + Cx^2 + Dx + E
        var solution:Array = new Array(4);

        // Using Ferrari's formula: http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution
        var Alpha:Number = ((-3 * (B * B)) / (8 * (A * A))) + (C / A);
        var Beta:Number = ((B * B * B) / (8 * A * A * A)) - ((B * C) / (2 * A * A)) + (D / A);          
        var Gamma:Number = ((-3 * B * B * B * B) / (256 * A * A * A * A)) + ((C * B * B) / (16 * A * A * A)) - ((B * D) / (4 * A * A)) + (E / A);

        var P:Number = ((-1 * Alpha * Alpha) / 12) - Gamma; 
        var Q:Number = ((-1 * Alpha * Alpha * Alpha) / 108) + ((Alpha * Gamma) / 3) - ((Beta * Beta) / 8);

        var PreRoot1:Number = ((Q * Q) / 4) + ((P * P * P) / 27);
        var R:ComplexNumber = ComplexNumber.add(new ComplexNumber((-1 * Q) / 2), ComplexNumber.sqrt(new ComplexNumber(PreRoot1)));

        var U:ComplexNumber = ComplexNumber.pow(R, 1/3);

        var preY1:Number = (-5 / 6) * Alpha;
        var RedundantY:ComplexNumber = ComplexNumber.add(new ComplexNumber(preY1), U);

        var Y:ComplexNumber;

        if(U.isZero()){
            var preY2:ComplexNumber = ComplexNumber.pow(new ComplexNumber(Q), 1/3);

            Y = ComplexNumber.subtract(RedundantY, preY2);
        } else{
            var preY3:ComplexNumber = ComplexNumber.multiply(new ComplexNumber(3), U);
            var preY4:ComplexNumber = ComplexNumber.divide(new ComplexNumber(P), preY3);

            Y = ComplexNumber.subtract(RedundantY, preY4);
        }

        var W:ComplexNumber = ComplexNumber.sqrt(ComplexNumber.add(new ComplexNumber(Alpha), ComplexNumber.multiply(new ComplexNumber(2), Y)));

        var Two:ComplexNumber = new ComplexNumber(2);
        var NegativeOne:ComplexNumber = new ComplexNumber(-1);

        var NegativeBOverFourA:ComplexNumber = new ComplexNumber((-1 * B) / (4 * A));
        var NegativeW:ComplexNumber = ComplexNumber.multiply(W, NegativeOne);

        var ThreeAlphaPlusTwoY:ComplexNumber = ComplexNumber.add(new ComplexNumber(3 * Alpha), ComplexNumber.multiply(new ComplexNumber(2), Y));
        var TwoBetaOverW:ComplexNumber = ComplexNumber.divide(new ComplexNumber(2 * Beta), W);

        solution["root1"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.add(W, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.add(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root2"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.subtract(NegativeW, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.subtract(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root3"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.subtract(W, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.add(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root4"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.add(NegativeW, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.subtract(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));

        return solution;
    }

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

Например, это уравнение: у = 0.9604000000000001x ^ 4 - 5.997600000000001x ^ 3 + 13.951750054511718x ^ 2 - 14.326264455924333x + 5,474214401412618

Возвращает корни: 1,7820304835380467 + 0i +1,34041662585388 + 0i 1.3404185025061823+ 0i 1.7820323472855648 + 0i

Если я построю график этого конкретного уравнения, я могу видеть, что фактические корни ближе к 1,2 и 2,9 (приблизительно).Я не могу отклонить четыре неправильных корня как случайные, потому что они на самом деле являются двумя корнями для первой производной уравнения:

y = 3.8416x ^ 3 - 17.9928x ^ 2 + 27.9035001x - 14.326264455924333

Имейте в виду, что я на самом деле не ищу конкретные корни уравнения, которое я опубликовал.У меня вопрос, есть ли какой-то особый случай, который я не принимаю во внимание.

Есть идеи?

Answer 1

Для нахождения корней многочленов степени> = 3 у меня всегда были лучшие результаты, используя Дженкинса-Трауба (http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm), чем явные формулы.

Answer 2

Я не знаю, почему решение Ferrari не работает, но я попытался использовать стандартный численный метод (создать сопутствующую матрицу и вычислить ее собственные значения), и я получил правильное решение, то есть два действительных корня в 1,2 и 1,9 .

Этот метод не для слабонервных. После построения вспомогательной матрицы полинома вы запускаете QR алгоритм , чтобы найти собственные значения этой матрицы. Это нули многочлена.

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

Answer 3

Другие ответы - хороший и здравый совет. Однако, вспоминая мой опыт внедрения метода Ferrari в Forth, я думаю, что ваши неправильные результаты, вероятно, вызваны: 1. неправильной реализацией необходимых и довольно хитрых комбинаций знаков, 2. еще не осознавая, что «.. == beta» с плавающей точкой должно стать "abs (.. - beta)

Для этой конкретной проблемы мой код Forth в диагностическом режиме возвращает:

x1 =  1.5612244897959360787072371026316680470492e+0000 -1.6542769593216835969789894020584464029664e-0001 i
 --> -4.2123274051525879873007970023884313331788e-0054  3.4544674220377778501545407451201598284464e-0077 i
x2 =  1.5612244897959360787072371026316680470492e+0000  1.6542769593216835969789894020584464029664e-0001 i
 --> -4.2123274051525879873007970023884313331788e-0054 -3.4544674220377778501545407451201598284464e-0077 i
x3 =  1.2078440724224197532447709413299479764843e+0000  0.0000000000000000000000000000000000000000e-0001 i
 --> -4.2123274051525879873010733597821943554068e-0054  0.0000000000000000000000000000000000000000e-0001 i
x4 =  1.9146049071693819497220585618954851525216e+0000 -0.0000000000000000000000000000000000000000e-0001 i
 --> -4.2123274051525879873013497171759573776348e-0054  0.0000000000000000000000000000000000000000e-0001 i

Текст после "->" следует из обратной подстановки корня в исходное уравнение.

Для справки, вот результаты Mathematica / Alpha с максимально возможной точностью, которые мне удалось установить:

Mathematica:
x1 = 1.20784407242
x2 = 1.91460490717
x3 = 1.56122449 - 0.16542770 i
x4 = 1.56122449 + 0.16542770 i 

Answer 4

Вы можете увидеть мой ответ на связанный вопрос . Я поддерживаю точку зрения Оливье: путь может быть просто подход сопутствующей матрицы / собственного значения (очень стабильный, простой, надежный и быстрый).

Редактировать

Полагаю, для меня не повредит, если я воспроизведу здесь ответ для удобства:

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

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

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

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

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

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

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

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

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

Удачи.

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

Answer 5

Хорошей альтернативой уже упомянутым методам является Алгоритм TOMS 326 , основанный на статье "Корни полиномов низкого порядка", подготовленной Теренсом RFNonweiler CACM (апрель 1968 г.).* Это алгебраическое решение полиномов 3-го и 4-го порядка, которое является достаточно компактным, быстрым и достаточно точным.Это намного проще, чем Дженкинс Трауб.

Имейте в виду, однако, что код TOMS работает не очень хорошо.

На этой странице Iowa Hills Root Solver есть код для поиска четвертого / кубического корня,немного более изысканным.У этого также есть искатель корня типа Дженкинса Трауба.

Answer 6

Я также реализовал эту функцию и получил вменяемый результат, как и вы.После некоторой отладки необходимо пройти проверку условия "\ beta <0", чего не произошло. </p>