一元三次方程的实数根

一元三次方程的实数根

一元三次方程的一般形式

一元三次方程的一般形式为

ax3+bx2+cx+d=0,a≠0.

令x=y−b3a，则原式变成

a(y−b3a)3+b(y−b3a)2+c(y−b3a)+d=0,a(y3−by2a+b2y3a2−b327a3)+b(y2−2by3a+b29a2)+c(y−b3a)+d=0,ay3−by2+b23ay−b327a2+by2−2b23ay+b39a2+cy−bc3a+d=0,ay3+(c−b23a)y+(d+2b327a2−bc3a)=0,y3+(ca−b23a2)y+(da+2b327a3−bc3a2)=0.

如此一来二次项就不见了，化成y3+3py+2q=0，其中p=c3a−b29a2，q=d2a+b327a3−bc6a2。

方程y3+3py+2q=0的实数根

令Δ=q2+p3，则

(i) 当Δ>0时，方程只有一个实根。

y=−q+Δ−−√−−−−−−−−√3+−q−Δ−−√−−−−−−−−√3.

(ii) 当Δ=0时，方程有三个实根，其中至少有两个相等的实根。

y1=−2q√3,y2=y3=q√3.

(iii) 当Δ<0时，方程有三个不等实根。

令α=13arccos−q−p−−−√p2，则

y1=2−p−−−√cosα,y2=2−p−−−√cos(α+120∘),y3=2−p−−−√cos(α+240∘).

C++代码实现

// π
const double PIE = 3.1415926535897932384626433832795;
// 求一个实数的立方根
static double SolveCubicRoot(double value);
// 求一元一次方程的实根：ax+b=0
static std::vector<double> SolveLinearEquation(double a, double b);
// 求一元二次方程的实根：ax^2+bx+c=0
static std::vector<double> SolveQuadraticEquation(double a, double b, double c);

// 求一元三次方程的实根：ax^3+bx^2+cx+d=0
std::vector<double> SolveCubicEquation(double a, double b, double c, double d)
{
    // 判断三次项系数是否为零
    if (a == 0)
    {
        return SolveQuadraticEquation(b, c, d);
    }

    std::vector<double> root;

    // 系数
    double p = (c / a - b * b / (3 * a * a)) / 3;
    double q = (d / a + 2 * b * b * b / (27 * a * a * a) - b * c / (3 * a * a)) / 2;
    double diff = -b / (3 * a);

    // 判别式
    double delta = p * p * p + q * q;
    if (delta > 0)
    {
        // 方程只有一个实根
        double sqrtDelta = sqrt(delta);
        root.push_back(SolveCubicRoot(-q + sqrtDelta) + SolveCubicRoot(-q - sqrtDelta));
    }
    else if (delta < 0)
    {
        // 方程有三个不等的实根
        double angle = acos(-q * sqrt(-p) / (p * p)) / 3;
        root.push_back(2.0 * sqrt(-p) * cos(angle));
        root.push_back(2.0 * sqrt(-p) * cos(angle + 2 * PIE / 3));
        root.push_back(2.0 * sqrt(-p) * cos(angle + 4 * PIE / 3));
    }
    else
    {
        // 方程有三个实根，其中至少有两个相等的实根
        if (q == 0)
        {
            root.push_back(0);
        }
        else
        {
            root.push_back(SolveCubicRoot(q));
            root.push_back(-2 * SolveCubicRoot(q));
        }
    }

    // 将结果加上余量
    for (int i = 0, maxSize = root.size(); i < maxSize; ++i)
    {
        root[i] += diff;
    }

    return root;
}

double SolveCubicRoot(double value)
{
    return value > 0 ? pow(value, 1.0 / 3) : -pow(-value, 1.0 / 3);
}

std::vector<double> SolveLinearEquation(double a, double b)
{
    std::vector<double> root;

    // 判断一次项系数是否为零
    if (a != 0)
    {
        root.push_back(-b / a);
    }

    return root;
}

std::vector<double> SolveQuadraticEquation(double a, double b, double c)
{
    // 判断二次项系数是否为零
    if (a == 0)
    {
        return SolveLinearEquation(b, c);
    }

    std::vector<double> root;

    // 计算判别式
    double delta = b * b - 4 * a * c;
    if (delta > 0)
    {
        // 方程有两个不等的实根
        root.push_back((-b + sqrt(delta)) / (2 * a));
        root.push_back((-b - sqrt(delta)) / (2 * a));
    }
    else if (delta == 0)
    {
        // 方程有两个相等的实根
        root.push_back(-b / (2 * a));
    }

    return root;
}