Specialised algorithm to find positive real solutions to quartic equations?

user_789308 0
Viewed 3853

I'm looking for a specialised algorithm to find positive real solutions to quartic equations with real coefficients (also know as bi-quadratic or polynomial equations of order 4). They have the form:

a4 x4 + a3 x3 +a2 x2 +a1 x + a0 = 0

with a1, a2,... being real numbers.

It's supposed to run on a microcontroller, which will need to do quite a lot of those calculations. So performance is an issue. That's why I'm looking for a specialised algorithm for positive solutions. If possible I'd like it to compute the exact solutions.

I know there is a general way to compute the solution of a quartic equation but it is rather involved in terms of computation.

Can someone point me in the right direction?

Edit:

Judging from the answers: Some seem to have misunderstood me (though I was pretty clear about it). I know of the standard ways of solving quartic equations. They don't do it for me - neither they fit in the memory nor are they sufficiently fast. What I would need is a high accuracy highly efficient algorithm to find only real solutions (if that helps) to quartic equations with real coefficients. I'm not sure there is such an algorithm, but I thought you guys might know. P.S.: The downvotes didn't come from me.

microcontroller math algorithm
7 Answers

Here is the C/C++ code that I wrote. But it only gives you real roots of the equation

Here is the C/C++ code that I wrote. But it only gives you real roots of the equation

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

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

        return nroots;
    }

    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;
}
user_21668570
Related