Специализированный алгоритм поиска положительных вещественных решений для квартичных уравнений?

Question

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

a ₄ x ⁴ + a ₃ x ³ + a ₂ x ² + a ₁ x + a ₀ = 0

с ₁ , a ₂ , ... являющиеся действительными числами.

Предполагается, что он будет работать на микроконтроллере, который должен будет выполнять довольно много таких вычислений.Так что производительность - это проблема.Вот почему я ищу специализированный алгоритм для положительных решений.Если возможно, я бы хотел, чтобы он вычислял точные решения.

Я знаю, что существует общий способ вычисления решения квартического уравнения, но он скорее связан с вычислениями.

Может кто-нибудь направить меня в правильном направлении?

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

Судя по ответам: некоторые, кажется, меня неправильно поняли (хотя я был совершенно уверен в этом). Мне известны стандартные способы решения квартичных уравнений.Они не делают это для меня - они не умещаются в памяти и не достаточно быстры.Что мне нужно, так это высокоточный высокоэффективный алгоритм, позволяющий находить только реальные решения (если это помогает) квартичных уравнений с реальными коэффициентами.Я не уверен, что есть такой алгоритм, но я думал, что вы, ребята, могли бы знать.PS: от меня не исходили голоса.

Answer 1

Это одна из тех ситуаций, когда, вероятно, легче найти все корни, используя сложную арифметику, чем только найти положительные реальные корни.И поскольку кажется, что вам нужно найти несколько корней одновременно, я бы рекомендовал использовать метод Дюранда-Кернера, который в основном является усовершенствованием метода Вейерштрасса:

http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method

Метод Вейерштрасса, в свою очередь, является усовершенствованным методом Ньютона, который решает параллельно для всех корней многочлена (и он имеет большое преимущество в том, что его просто убить мозгами).В целом он сходится примерно с квадратичной скоростью, хотя только для нескольких корней.Для большинства квартичных полиномов вы можете в значительной степени прибить корни всего за несколько итераций.Если вам нужно решение более общего назначения, вам следует использовать вместо него Jenkins-Traub:

http://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_method

Это быстрее для полиномов более высокой степени и в основном работает путем преобразования задачи вНахождение собственных значений сопутствующей матрицы:

http://en.wikipedia.org/wiki/Companion_matrix

РЕДАКТИРОВАТЬ: В качестве второго предложения, вы также можете попробовать использовать метод мощности на сопутствующей матрице.Поскольку ваше уравнение имеет только неотрицательные коэффициенты, может оказаться полезным применить теорему Перрона-Фробениуса к сопутствующей матрице.Как минимум, это свидетельствует о том, что существует хотя бы один неотрицательный корень:

http://en.wikipedia.org/wiki/Perron%E2%80%93Frobenius_theorem

Answer 2

Да, есть общие способы.Вам нужен алгоритм поиска корня, такой как брекетинг и деление на две части, секущий, ложная позиция, Риддер, Ньютон-Рафсон, дефляция, Мюллер, Лагерр или Дженкинс-Трауб - я кого-нибудь пропустил?

Проверить "ЧисловойРецепты "для деталей.

Answer 3

Вот код C / C ++, который я написал. Но это только дает вам реальные корни уравнения

#include <stdio.h>
#include <iostream>
#include <math.h>
/*--------------------------------------------

 --------------------------------------------*/
double cubic(double b,double c,double d)
{
    double p=c-b*b/3.0;
    double q=2.0*b*b*b/27.0-b*c/3.0+d;

    if(p==0.0) return pow(q,1.0/3.0);
    if(q==0.0) return 0.0;

    double t=sqrt(fabs(p)/3.0);
    double g=1.5*q/(p*t);
    if(p>0.0)
    return -2.0*t*sinh(asinh(g)/3.0)-b/3.0;


    if(4.0*p*p*p+27.0*q*q<0.0)
    return 2.0*t*cos(acos(g)/3.0)-b/3.0;

    if(q>0.0)
    return -2.0*t*cosh(acosh(-g)/3.0)-b/3.0;

    return 2.0*t*cosh(acosh(g)/3.0)-b/3.0;
}
/*--------------------------------------------

 --------------------------------------------*/
int quartic(double b,double c,double d,double e,double* ans)
{

    double p=c-0.375*b*b;
    double q=0.125*b*b*b-0.5*b*c+d;
    if(q==0)
    {
        double m=0.25*p*p+0.01171875*b*b*b*b-e+0.25*b*d-0.0625*b*b*c;
        if(m<0.0) return 0;
        int nroots=0;

        if(m-0.5*p>=0.0)
        {
            double absy=sqrt(m-0.5*p);
            ans[nroots++]=absy-0.25*b;
            ans[nroots++]=-absy-0.25*b;
        }
        if(m+0.5*p>=0.0)
        {
            double absy=sqrt(m+0.5*p);
            ans[nroots++]=absy-0.25*b;
            ans[nroots++]=-absy-0.25*b;
        }

        return nroots;
    }
    double m=cubic(p,0.25*p*p+0.01171875*b*b*b*b-e+0.25*b*d-0.0625*b*b*c,-0.125*q*q);

    if(m<0.0) return 0;
    double sqrt_2m=sqrt(2.0*m);
    int nroots=0;
    if(-m-p+q/sqrt_2m>=0.0)
    {
        double delta=sqrt(2.0*(-m-p+q/sqrt_2m));
        ans[nroots++]=0.5*(-sqrt_2m+delta)-0.25*b;
        ans[nroots++]=0.5*(-sqrt_2m-delta)-0.25*b;
    }

    if(-m-p-q/sqrt_2m>=0.0)
    {
        double delta=sqrt(2.0*(-m-p-q/sqrt_2m));
        ans[nroots++]=0.5*(sqrt_2m+delta)-0.25*b;
        ans[nroots++]=0.5*(sqrt_2m-delta)-0.25*b;
    }

    return nroots;
}
/*--------------------------------------------

 --------------------------------------------*/
int main(int nargs,char* args[])
{
    if(nargs!=6)
    {
        printf("5 arguments are needed\n");
        return EXIT_FAILURE;
    }
    double a=atof(args[1]);
    double b=atof(args[2]);
    double c=atof(args[3]);
    double d=atof(args[4]);
    double e=atof(args[5]);
    if(a==0.0)
    {
        printf("1st argument should be nonzero\n");
        return EXIT_FAILURE;
    }

    int nroots;
    double ans[4];
    nroots=quartic(b/a,c/a,d/a,e/a,ans);
    if(nroots==0)
        printf("Equation has no real roots!\n");
    else
    {
        printf("Equation has %d real roots: ",nroots);
        for(int i=0;i<nroots-1;i++) printf("%.16lf, ",ans[i]);
        printf("%.16lf\n",ans[nroots-1]);
    }

    return EXIT_SUCCESS;
}

Answer 4

Я понимаю, что этот ответ довольно поздно, но я думаю, что хорошей альтернативой уже упомянутым методам является Алгоритм TOMS 326 , основанный на статье Теренса Р.Ф. "Корни полиномов низкого порядка" Nonweiler CACM (апрель 1968).

Это алгебраическое решение полиномов 3-го и 4-го порядка, которое является достаточно компактным и быстрым. Это, конечно, намного проще и намного быстрее, чем Jenkins Traub, и не требует итерации.

Однако не используйте код TOMS, так как он был сделан довольно плохо.

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

Answer 5

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

Вам повезло, что аналитическое решение вообще есть.Если бы вы были полиномом 5-го порядка, этого бы даже не было.

Answer 6

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

Answer 7

Можете ли вы предоставить хорошие начальные значения, чтобы всегда находить все решения.Метод Ньютона быстро сходится.

Я проверил в Maxima:

solve(a*x^4+b*x^3+c*x^2+d*x+c=0,x);

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

Однако, если коэффициенты постоянны, вы можете просто реализовать прямую формулу.Вы можете получить решение, установив максимумы или введя уравнение на wolframalpha.com