求解一元高次多项式方程的所有实数根的算法

一元高次多项式的形式可以表示为：
$$f(x)=c_0+c_{1}x+c_2{x}^2+c_3{x}^3+...+c_n{x}^n=0$$

一般来说，一元高次多项式方程很难有显示解析解，大多使用数值算法计算方程的根，例如牛顿法、割线法等。但这些算法的缺点是计算得到的解受到迭代初始值的影响，另外在某些特定情况下，迭代也有可能不收敛。对于拥有多个实数根的多项式方程，使用这些迭代算法难以获得方程所有实数根。
本文将介绍一种新的算法，通过递归降次，将高次多项式方程的求解问题逐步转化为可以解析求解的一元二次方程，进而求解出多项式的所有实数根。

原理

首先，为了求解多项式方程的所有实根，我们需要对多项式函数进行分段。使用$x_1,x_2,x_3,...,x_k$这一系列的点将实数轴分成若干区间段，以使得在每一个区间段上函数f(x)都是单调的。
当我们获得了这一系列的分隔点后，例如上述k个分隔点，将实数轴分为k+1个区间段，（包括$(-\infty ,x_1]$和$[x_k,+\infty )$两个区间）。在每一个区间中函数都是单调的。对于区间$[x_i,x_{i+1}]$，如果$f(x_i)$和$f(x_{i+1})$符号相反，由于函数在区间中是单调的，在该区间中必然有且仅有一个实根。这时只要我们在该区间中使用最简单的二分法就能得到该区间中的或者一个实根。
当然，如果$f(x_i)$和$f(x_{i+1})$符号相同，由于函数单调，在该区间中必然不存在方程的实根。
对于$(-\infty ,x_1]$这个区间，可以先算出$f(x_1)$和$f(x_1-1)$，如果$f(x_1-1)$和$f(x_1)$同号，且$f(x_1-1)$比$f(x_1)$更远离0，由于函数的单调性，在这个区间上必然不存在实根。如果不是这种情况，我们就需要依次计算$f(x_1-1),f(x_1-2),f(x_1-4),...,f(x_1-2^m),...$，直到找到一个与$f(x_1)$符号相反的值，例如$x_1-2^m$，然后再在区间$[x_1-2^m,x_1]$上使用二分法求解出相应的实数根。
对于区间$[x_k,+\infty )$也使用一样的处理方法，只是将上面的减号变成加号即可。
使用以上的方法，我们发现，只要我们能够找到这一系列的分隔点，使得分隔出的每一个区间上多项式函数都是单调的，那么这个多项式方程求解的问题就解决了。
那么问题就转换为如何找到这一系列的分隔点。对于多项式函数，函数的一阶导数等于0的点，正好就是这样的分隔点。我们只要求解出函数f(x)的所有一阶导数等于0的点即可。即问题转化为求解方程$$f^{{}^{\prime }}(x)=0$$也就是
$$f^{{}^{\prime }}(x)=c_1+2c_{2}x+3c_3{x}^2+4c_4{x}^3+...+n{c}_n{x}^{n-1}=0$$

经过以上的处理，我们将求解f(x)=0方程的问题转换成了求解其一阶导函数等于0的方程求解问题，这样使得我们求解的方程的次数下降了一阶。
使用这种方法，递归处理，我们最终一定能够将问题逐步转换为可以解析求解的一元二次方程。
可以看到，以上求解的过程完全不需要外部指定任何初始值以及其他变量，只需要输入方程每一项的系数，就可以求解出所有了方程实根。

代码实现（C#）

    /// <summary>
    /// 求解单一未知数多项式方程，c0+c1*x^1+c2*x^2+c3*x^3+...+cn*x^n=0
    /// </summary>
    public class SingleUnknownPolynomialSolver
    {
        /// <summary>
        /// 传入系数数组，构建实例
        /// </summary>
        /// <param name="coefs">方程系数：c0、c1、c2、c3、...、cn</param>
        public SingleUnknownPolynomialSolver(params double[] coefs)
        {
            if (coefs == null) _c = null;
            int cLength = 0;
            for (int i = 0; i < coefs.Length; i++)
            {
                if (coefs[coefs.Length - 1 - i] == 0) continue;
                cLength = coefs.Length - i;
                break;
            }
            _c = coefs.Take(cLength).ToArray();
        }


        private double[] _c;
        private Mathematics.ApproxEqual _approx = new Mathematics.ApproxEqual();

        /// <summary>
        /// 获取或设置近似相等比较器，近似相等比较器不可以为null，如果设置null则使用默认的近似相等比较器
        /// </summary>
        public Mathematics.ApproxEqual ApproxEqual
        {
            get => _approx;
            set
            {
                if (value == null) _approx = new Mathematics.ApproxEqual();
                else _approx = value;
            }
        }

        /// <summary>
        /// 求解方程的所有实数解
        /// </summary>
        /// <returns>返回所有实数解，无解返回空数组，无穷解返回null</returns>
        public double[] Solve()
        {
            if (_c == null || _c.Length == 0) return null;
            if (_c.Length == 1) return new double[] { };
            if (_c.Length == 2) return new double[] { -_c[0] / _c[1] };
            if (_c.Length == 3)
            {
                double delta = _c[1] * _c[1] - 4 * _c[0] * _c[2];
                if (delta < 0) return new double[] { };
                double sqrtDelta = Math.Sqrt(delta);
                if (_approx.IsEqual(sqrtDelta / _c[2], 0)) return new double[] { -_c[1] / _c[2] / 2 };
                return new double[] { (-_c[1] + sqrtDelta) / _c[2] / 2, (-_c[1] - sqrtDelta) / _c[2] / 2 };
            }
            var partialSolver = PartialSolver();
            var partialSolutions = partialSolver.Solve();
            if (partialSolutions.Length == 0)
            {
                double? x_down = OppositePoint(0, false);
                if (x_down != null)
                {
                    double? solution = SolvePoint(x_down ?? 0, 0);
                    if (solution != null) return new double[] { solution ?? 0 };
                    else return new double[] { };
                }
                double? x_up = OppositePoint(0, true);
                if (x_up != null)
                {
                    double? solution = SolvePoint(0, x_up ?? 0);
                    if (solution != null) return new double[] { solution ?? 0 };
                    else return new double[] { };
                }
                return new double[] { };
            }
            List<double> solutionList = new List<double>();
            var partialSolutionList = partialSolutions.ToList();
            partialSolutionList.Sort();
            partialSolutions = partialSolutionList.ToArray();
            double? x_d = OppositePoint(partialSolutions.First(), false);
            if (x_d != null)
            {
                double? solution = SolvePoint(x_d ?? 0, partialSolutions.First());
                if (solution != null) solutionList.Add(solution ?? 0);
            }
            double? x_u = OppositePoint(partialSolutions.Last(), true);
            if (x_u != null)
            {
                double? solution = SolvePoint(partialSolutions.Last(), x_u ?? 0);
                if (solution != null) solutionList.Add(solution ?? 0);
            }
            if (partialSolutions.Length >= 2)
            {
                for (int i = 1; i < partialSolutions.Length; i++)
                {
                    double? solution = SolvePoint(partialSolutions[i - 1], partialSolutions[i]);
                    if (solution != null) solutionList.Add(solution ?? 0);
                }
            }
            solutionList.Sort();
            for (int i = 1; i < solutionList.Count; i++)
            {
                if (_approx.IsEqual(solutionList[i - 1], solutionList[i]))
                {
                    solutionList.RemoveAt(i);
                    i--;
                }
            }
            return solutionList.ToArray();
        }

        /// <summary>
        /// 计算函数值
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private double CalcValue(double x)
        {
            double y = 0;
            for (int i = 0; i < _c.Length; i++)
            {
                y += _c[i] * Math.Pow(x, i);
            }
            return y;
        }

        /// <summary>
        /// 一阶导函数方程
        /// </summary>
        /// <returns></returns>
        private SingleUnknownPolynomialSolver PartialSolver()
        {
            double[] coefs = new double[_c.Length - 1];
            for (int i = 1; i < _c.Length; i++)
            {
                coefs[i - 1] = _c[i] * i;
            }
            SingleUnknownPolynomialSolver solver = new SingleUnknownPolynomialSolver(coefs);
            solver.ApproxEqual = _approx;
            return solver;
        }

        /// <summary>
        /// 求解x点函数值符号相反的点，要求指定的方向上函数是单调的
        /// </summary>
        /// <param name="x"></param>
        /// <param name="upper">向上搜索还是向下搜索</param>
        /// <returns></returns>
        private double? OppositePoint(double x, bool upper = false)
        {
            double s = upper ? 1 : -1;
            double y = CalcValue(x);
            if (y == 0) return x;
            double m = 1;
            while (m < double.MaxValue)
            {
                double xp = x + m * s;
                double yp = CalcValue(xp);
                if (yp == 0) return xp;
                if (yp == y) return null;
                if (yp > y && y > 0) return null;
                if (yp < y && y < 0) return null;
                if (yp > 0 && y < 0) return xp;
                if (xp < 0 && y > 0) return xp;
                m = m * 2;
            }
            return null;
        }

        /// <summary>
        /// 求解两个点中间的方程解，要求在两点之间函数单调
        /// </summary>
        /// <param name="xl">下方点</param>
        /// <param name="xu">上方点</param>
        /// <returns></returns>
        private double? SolvePoint(double xl, double xu)
        {
            if (xl > xu) return null;
            if (xl == xu)
            {
                if (CalcValue(xl) == 0) return xl;
                else return null;
            }
            double l = xl;
            double u = xu;
            while (!_approx.IsEqual(l, u))
            {
                double yl = CalcValue(l);
                if (yl == 0) return l;
                double yu = CalcValue(u);
                if (yu == 0) return u;
                if (yl > 0 && yu > 0) return null;
                if (yl < 0 && yu < 0) return null;
                double m = (l + u) / 2;
                double ym = CalcValue(m);
                if (ym == 0) return m;
                if (ym > 0 && yl < 0) u = m;
                else if (ym > 0 && yl > 0) l = m;
                else if (ym < 0 && yl < 0) l = m;
                else u = m;
            }
            return (u + l) / 2;
        }

    }

关于代码的一些说明

1、算法的核心在于Solve方法中，采用嵌套调用的形式。OppositePoint方法用于处理边界区间搜索异号函数值点，SolvePoint方法用于求解在指定的区间上使用二分法求出区间上的方程根。
2、在构造函数中对输入的系数数组进行了预处理，以保证多项式的最高次项的系数不为0。
3、上述代码中没有给出Mathematics.ApproxEqual类的具体定义。实际上ApproxEqual是一个近似相等比较器，它的IsEqual方法可以比较两个实数在指定的精度和算法下是否近似相等。在实际的应用中，读者可以根据自己的需要定义自己的近似相等比较器的具体算法和参数。
4、完整的代码可以下载：https://download.csdn.net/download/weixin_45107579/13712341