0.引言

在学习过程中需要求解一个一元四次方程。
由于带有未知参数，本想使用MATLAB求解出它的解析解然后写入C++进行计算。后来发现求解出来的解太长。

要是直接放入C++中，通过解析解来实时求它的解甚至会比求解一元四次方程花费的时间更多。于是另谋出路，查阅资料如何求解一元四次方程。

1.方法一

Ferrari方法，这也是最容易Google到的方法。

• wiki

时间紧张原理就不看了，直接找实现。巧了，某百科上就有code.Link.

Ferrari.cpp

#include <math.h>
#include <float.h>
#include <complex>
#include <iostream>
/******************************************************************************\
对一个复数 x 开 n 次方
\******************************************************************************/
std::complex<double>  sqrtn(const std::complex<double>&x,double n)
{
    double r = hypot(x.real(),x.imag()); //模
    if(r > 0.0)
    {
        double a = atan2(x.imag(),x.real()); //辐角
        n = 1.0 / n;
        r = pow(r,n);
        a *= n;
        return std::complex<double>(r * cos(a),r * sin(a));
    }
    return std::complex<double>();
}
/******************************************************************************\
使用费拉里法求解一元四次方程 a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
\******************************************************************************/
std::complex<double> Ferrari(std::complex<double> x[4]
,std::complex<double> a
,std::complex<double> b
,std::complex<double> c
,std::complex<double> d
,std::complex<double> e)
{
    a = 1.0 / a;
    b *= a;
    c *= a;
    d *= a;
    e *= a;
    std::complex<double> P = (c * c + 12.0 * e - 3.0 * b * d) / 9.0;
    std::complex<double> Q = (27.0 * d * d + 2.0 * c * c * c + 27.0 * b * b * e - 72.0 * c * e - 9.0 * b * c * d) / 54.0;
    std::complex<double> D = sqrtn(Q * Q - P * P * P,2.0);
    std::complex<double> u = Q + D;
    std::complex<double> v = Q - D;
    if(v.real() * v.real() + v.imag() * v.imag() > u.real() * u.real() + u.imag() * u.imag())
    {
        u = sqrtn(v,3.0);
    }
    else
    {
        u = sqrtn(u,3.0);
    }
    std::complex<double> y;
    if(u.real() * u.real() + u.imag() * u.imag() > 0.0)
    {
        v = P / u;
        std::complex<double> o1(-0.5,+0.86602540378443864676372317075294);
        std::complex<double> o2(-0.5,-0.86602540378443864676372317075294);
        std::complex<double>&yMax = x[0];
        double m2 = 0.0;
        double m2Max = 0.0;
        int iMax = -1;
        for(int i = 0;i < 3;++i)
        {
            y = u + v + c / 3.0;
            u *= o1;
            v *= o2;
            a = b * b + 4.0 * (y - c);
            m2 = a.real() * a.real() + a.imag() * a.imag();
            if(0 == i || m2Max < m2)
            {
                m2Max = m2;
                yMax = y;
                iMax = i;
            }
        }
        y = yMax;
    }
    else
    {//一元三次方程，三重根
        y = c / 3.0;
    }
    std::complex<double> m = sqrtn(b * b + 4.0 * (y - c),2.0);
    if(m.real() * m.real() + m.imag() * m.imag() >= DBL_MIN)
    {
        std::complex<double> n = (b * y - 2.0 * d) / m;
        a = sqrtn((b + m) * (b + m) - 8.0 * (y + n),2.0);
        x[0] = (-(b + m) + a) / 4.0;
        x[1] = (-(b + m) - a) / 4.0;
        a = sqrtn((b - m) * (b - m) - 8.0 * (y - n),2.0);
        x[2] = (-(b - m) + a) / 4.0;
        x[3] = (-(b - m) - a) / 4.0;
    }
    else
    {
        a = sqrtn(b * b - 8.0 * y,2.0);
        x[0] =
        x[1] = (-b + a) / 4.0;
        x[2] =
        x[3] = (-b - a) / 4.0;
    }
return x[4];
}

int main()
{
std::complex<double> x[4];
x[4] = Ferrari(x,1,2,3,4,5); //验证费拉里法
std::cout<<"root1: "<<x[0]<<std::endl<<"root2: "<<x[1]<<std::endl<<"root3: "<<x[2]<<std::endl<<"root4: "<<x[3]<<std::endl;
return true;
}

测试系数：1,2,3,4,5

确保结果正确，使用matlab验证一下：

p = [1 2 3 4 5];
x = roots(p)

结果一致，OK。

2.方法二

使用多项式的伴随矩阵进行求解。

• ref

对于多项式：$P(x)=x^n+c_{n-1}{x}^{n-1}+c_{n-2}{x}^{n-2}+\cdots +c_{1}x+c_0$
其伴随矩阵：

$$M_x=\left[\begin{array}{ccccc}0& 0& ...& 0& -c_0\\ 1& 0& ...& 0& -c_1\\ 0& 1& ...& 0& -c_2\\ ...& ...& ...& ...& ...\\ 0& 0& ...& 1& -c_{n-1}\end{array}\right]$$
的特征值就是$P(x)$的根。

于是使用Eigen进行多项式的求解：

#include <iostream>
#include <Eigen/Core>
// 稠密矩阵的代数运算（逆，特征值等）
#include <Eigen/Dense>
using namespace std;

int main( int argc, char** argv )
{

    Eigen::Matrix<double, 4, 4> matrix_44;
    //复数动态矩阵
	Eigen::Matrix<complex<double>, Eigen::Dynamic, Eigen::Dynamic> matrix_eigenvalues;
	//同样测试 12345
    matrix_44 << 0, 0, 0, -5,
				 1, 0, 0, -4,
				 0, 1, 0, -3,
				 0, 0, 1, -2;

	std::cout<<"matrix_44: "<<std::endl<<matrix_44<<std::endl<<std::endl;

	matrix_eigenvalues = matrix_44.eigenvalues();
	//std::cout<<"matrix_44.eigenvalues: "<<std::endl<<matrix_44.eigenvalues()<<std::endl;
	std::cout<<"matrix_eigenvalues: "<<std::endl<<matrix_eigenvalues<<std::endl;
	return 0;
}

结果与上面同样一致：

3.结论

数学功底越深厚，代码复杂度越低！！

• 一元二次方程

% 求解 ax2+bx+c=0   2次方程有两个个解

syms a b c x;
[x] = solve(a*x^2+b*x+c==0,x,'Real',false)
syms x
[x] = solve(x^2+x+1==0,x,'Real',false)  %false 求复数解 true 求实数解
%图形
x = -2:0.01:1;
y = x.^2+x+1;
y1 = zeros(size(x));
plot(x,y,'b');
hold on 
plot(x,y1,'r--')

• 一元三次方程

% 求解 ax^3+bx^2+cx+d=0
syms a b c d x
[x] = solve(a*x^3+b*x^2+c*x+d==0,x) 
%  root(az3+bz2+cz+d,z,1)
syms x;
[x] = solve(x^3+x^2+x+1==0,x)  % 解x = -1,-i,i
syms x;
[x] = solve(x^3+x^2+x+1==0,x,'real',true) % x = -1
syms x;
[x] = solve(2*x^3+5*x^2+x+1==3,x) % x = -1

%图形
x = -4:0.01:3;
y = 2*x.^3+5*x.^2+x+1;
y1 = zeros(size(x));
plot(x,y,'b');
hold on 
plot(x,y1,'r--')

x = 

• 一元四次方程

% 求解 ax^4+bx^3+cx^2+dx+e=0
syms a b c d e x;
[x] = solve(a*x^4+b*x^3+c*x^2+d*x+e==0,x)
syms  x;
[x] = solve(x^4+x^3+x^2+x+1==1,x,'Real',true) %实数解

%图形
x = -8:0.01:4;
y = 2*x.^4+17*x.^3+7*x.^2+8*x+1;
y1 = zeros(size(x));
plot(x,y,'b');
hold on 
plot(x,y1,'r--')

x =

x =

• 一元五次方程

% 求解 ax^5+bx^4+cx^3+dx^2+ex+f==0
syms a b c d e f x
[x] = solve(a*x^5+b*x^4+c*x^3+d*x^2+e*x+f==0,x)
syms  x
[x] = solve(x^5+x^4+x^3+x^2+x+1==0,x,'Real',true)

%图形
x = -6:0.01:4;
y = 2*x.^5+12*x.^4+8*x.^3+20*x.^2+42*x+1;
y1 = zeros(size(x));
plot(x,y,'b');
hold on 
plot(x,y1,'r--')

x =

x =

• 一元六次方程

% 求解 ax^6+bx^5+cx^4+dx^3+ex^2+fx+g=0
syms a b c d e f g x
[x] = solve(a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g==0,x)
syms x
[x] = solve(x^6+x^5+x^4+x^3+x^2+x+1==1,x,'Real',true)

%图形
x = -30:0.01:20;
y = x.^6+25*x.^5+20*x.^4+30*x.^3+40*x.^2+50*x+1;
y1 = zeros(size(x));
plot(x,y,'b');
hold on 
plot(x,y1,'r--')

x =

x =