C# 的多项式方程求解的基本实现

raymoodlmd

已于 2024-04-11 15:16:12 修改

阅读量972

点赞数 21

文章标签：算法 c# 数学建模 matlab

于 2024-03-29 16:35:20 首次发布

本文链接：https://blog.csdn.net/raymoodlmd/article/details/137143400

版权

已知：已有很多数据，根据生成曲线方程并求得y=0时其多个方程解。

例如：x = [-5,- 4,- 3,- 2,- 1,0,1,2]

y = [-4697,- 1374,- 251, - 8, 3, - 2, 1, 84]

1、求其最小二乘法 k 阶多项式的拟合方程。

2、根据拟合方程求其与横坐标交点。

步骤：1、我们要写一个根据已知数据，求拟合方程参数的方法。

2、已知方程后，要写一个根据方程求解的方法

先已知的方法有很多，例如穷举法、二分法、牛顿-莱布利兹法、梯度法、碰撞体积法等等，各有优劣。博主这次根据自己的使用场景使用的是牛顿法。

a、牛顿法比较重要的是使用前提条件，连续可导（博主的使用场景是必定满足的，故就不考虑了）（基本知识，不记得可以去百度一下）

b、牛顿法的初值判断（这个其实看似容易，如果想准确，其实并不简单）

c、求导算法（基本知识可以去百度一下）

d、根据求出的解带入原方程，以此判断是否需要继续迭代。

纪要

工具使用：VS2022，百度在线图形计算器Desmos

使用相关库：MathNet

一、获得拟合参数方法

        /// <summary>
        /// 将点 (x,y) 拟合为 k 阶多项式 y 的最小二乘法，返回其最佳拟合参数为[p0, p1, p2,..., pk] 数组 p0、p1 分别指0次系数、1次系数
        /// </summary>
        /// <param name="xdata">x数组</param>
        /// <param name="ydata">y数组</param>
        /// <param name="order">拟合次数/多项式阶数</param>
        /// <returns>拟合参数数组</returns>
        public static double[] PolynomialPara(double[] xdata, double[] ydata, int order)
        {
            return Fit.Polynomial(xdata, ydata, order);
        }

运行

static double[] Xdata = { -5, -4, -3, -2, -1, 0, 1, 2 };
static double[] Ydata = { -4697, -1374, -251, -8, 3, -2, 1, 84 };


double[] _para = PolynomialPara(Xdata,Ydata,5);

结果

根据拟合参数，其实对应的k阶多项式方程就显而易见了

二、根据自变量x值,求结果的方法

/// <summary>
/// 返回 自变量带入多项式方程的结果
/// </summary>
/// <param name="polynomialPara">多项式参数</param>
/// <param name="x">自变量x</param>
/// <returns>多项式值</returns>
public static double PolynomialResult(double[] polynomialPara, double x)
{
    double result = 0.0;
    for (int i = 0; i < polynomialPara.Length; i++)
    {
        result += (polynomialPara[i] * Math.Pow(x, i));
    }

    return result;
}

三、根据拟合方程参数，进行一阶求导

/// <summary>
/// 多项式的1阶求导
/// </summary>
/// <param name="polynomialPara">多项式的参数</param>
/// <returns></returns>
public static double[]? FirstDerivative(double[] polynomialPara)
{
    List<double>? FirstDerivativepara = new List<double>();
    for (int i = 0; i < polynomialPara.Length; i++)
    {
        FirstDerivativepara.Add(polynomialPara[i] * i);
        FirstDerivativepara.Remove(0);
    }

    return FirstDerivativepara.ToArray();
}

四、牛顿法初值获得

基本逻辑，当我们设置一个初值进行迭代的时候，最终会获得与该初值距离最近的解；意味着如果方程有多解，我们就必须在每个解的周围都需要至少1个设置的初值。

博主有使用过根据区间产生随机值来生成初值，但是由于随机值不可控，故放弃了这一设想；采用区间二分法来获得初值，分n次，产生2^n + 1数量的初值。

/// <summary>
/// 二分法生成数 ex: FirstDichotomy([0,10], 2) => [0, 3.5, 5, 7.5, 10]
/// </summary>
/// <param name="interval">区间</param>
/// <param name="maxtimes">二分法次数</param>
/// <returns>二分法所的数的数组</returns>
public static double[] FirstDichotomy(double[] interval, int maxtimes = 1)
{
    List<double> reslut = new List<double>();

    double min = interval.Min();
    double max = interval.Max();
    double stepvalue = (max - min) / Math.Pow(2, maxtimes);
    if (stepvalue > 0)
    {
        double locValue = min;
        while (locValue <= max)
        {
            reslut.Add(locValue);
            locValue += stepvalue;
        }
        return reslut.ToArray();
    }
    else
    {
        return reslut.ToArray();
    }

}

五、牛顿法根据初值与多项式参数，求解

/// <summary>
/// 牛顿法求方程解
/// </summary>
/// <param name="x0">初始值</param>
/// <param name="para">多项式参数</param>
/// <param name="times">当前迭代次数</param>
/// <param name="maxTimes">最大迭代次数</param>
/// <param name="accuracy">精度</param>
/// <returns></returns>
public static double NewtonFunc(double x0, double[] para, int maxTimes, int times = 1, double accuracy = 0.0001)
{
    double x1 = 0;

    x1 = x0 - PolynomialResult(para, x0) / PolynomialResult(FirstDerivative(para), x0);


    if (Math.Abs(PolynomialResult(para, x1)) < accuracy || times > maxTimes)
    {
        //Trace.WriteLine($"times= {times} x = {x1} y = {PolynomialResult(para, x1)}");
        return x1;
    }
    times++;
    return NewtonFunc(x1, para, maxTimes, times);
}

最终完整代码

using MathNet.Numerics;
using System.Diagnostics;

namespace Poly
{
    internal class Program
    {

        static double[] Xdata = { -5, -4, -3, -2, -1, 0, 1, 2 };
        static double[] Ydata = { -4697, -1374, -251, -8, 3, -2, 1, 84 };
        static double[] interval = { -100, 100 }; // 生成二分法初值的区间

        static void Main()
        {

            double[] _para = PolynomialPara(Xdata, Ydata, 5);// 生成拟合参数
            double[] firstx0 = FirstDichotomy(interval, maxtimes: 15);// 生成2^15 + 1 个interval范围内的二分初值

            double[] resultx = new double[firstx0.Length];

            for (int i = 0; i < firstx0.Length; i++)
            {

                var re = NewtonFunc(firstx0[i], _para, maxTimes: 40);
                resultx[i] = Math.Round(re, 3);
                //resultx[i] = re;
            }
            var reslut = resultx.Distinct().ToArray();
            Print(reslut);

        }





    /// <summary>
    /// 将点 (x,y) 拟合为 k 阶多项式 y 的最小二乘法，返回其最佳拟合参数为[p0, p1, p2,..., pk] 数组 p0、p1 分别指0次系数、1次系数
    /// </summary>
    /// <param name="xdata">x数组</param>
    /// <param name="ydata">y数组</param>
    /// <param name="order">拟合次数/多项式阶数</param>
    /// <returns>拟合参数数组</returns>
    static double[] PolynomialPara(double[] xdata, double[] ydata, int order)
        {
            return Fit.Polynomial(xdata, ydata, order);
        }

        /// <summary>
        /// 返回 自变量带入多项式方程的结果
        /// </summary>
        /// <param name="polynomialPara">多项式参数</param>
        /// <param name="x">自变量x</param>
        /// <returns>多项式值</returns>
        public static double PolynomialResult(double[] polynomialPara, double x)
        {
            double result = 0.0;
            for (int i = 0; i < polynomialPara.Length; i++)
            {
                result += (polynomialPara[i] * Math.Pow(x, i));
            }

            return result;
        }


        /// <summary>
        /// 多项式的1阶求导
        /// </summary>
        /// <param name="polynomialPara">多项式的参数</param>
        /// <returns></returns>
        public static double[]? FirstDerivative(double[] polynomialPara)
        {
            List<double>? FirstDerivativepara = new List<double>();
            for (int i = 0; i < polynomialPara.Length; i++)
            {
                FirstDerivativepara.Add(polynomialPara[i] * i);
                FirstDerivativepara.Remove(0);
            }

            return FirstDerivativepara.ToArray();
        }

        /// <summary>
        /// 二分法生成数 ex: FirstDichotomy([0,10], 2) => [0, 3.5, 5, 7.5, 10]
        /// </summary>
        /// <param name="interval">区间</param>
        /// <param name="maxtimes">二分法次数</param>
        /// <returns>二分法所的数的数组</returns>
        public static double[] FirstDichotomy(double[] interval, int maxtimes = 1)
        {
            List<double> reslut = new List<double>();

            double min = interval.Min();
            double max = interval.Max();
            double stepvalue = (max - min) / Math.Pow(2, maxtimes);
            if (stepvalue > 0)
            {
                double locValue = min;
                while (locValue <= max)
                {
                    reslut.Add(locValue);
                    locValue += stepvalue;
                }
                return reslut.ToArray();
            }
            else
            {
                return reslut.ToArray();
            }

        }

        /// <summary>
        /// 牛顿法求方程解
        /// </summary>
        /// <param name="x0">初始值</param>
        /// <param name="para">多项式参数</param>
        /// <param name="times">当前迭代次数</param>
        /// <param name="maxTimes">最大迭代次数</param>
        /// <param name="accuracy">精度</param>
        /// <returns></returns>
        public static double NewtonFunc(double x0, double[] para, int maxTimes, int times = 1, double accuracy = 0.0001)
        {
            double x1 = 0;

            x1 = x0 - PolynomialResult(para, x0) / PolynomialResult(FirstDerivative(para), x0);


            if (Math.Abs(PolynomialResult(para, x1)) < accuracy || times > maxTimes)
            {
                //Trace.WriteLine($"times= {times} x = {x1} y = {PolynomialResult(para, x1)}");
                return x1;
            }
            times++;
            return NewtonFunc(x1, para, maxTimes, times);
        }




        /// <summary>
        /// 提供打印列表工具
        /// </summary>
        /// <typeparam name="T">变量类型</typeparam>
        /// <param name="list"List 列表</param>
        static void Print<T>(IEnumerable<T> list)
        {
            Trace.Write("[");
            foreach (var i in list)
            {
                Trace.Write($"{i.ToString()},");

            }
            Trace.WriteLine("]");

        }
    }
}

运行结果

该结果比对图形计算器的可视结果，基本一致。

结论：

优点：利用了牛顿法收敛快的优点，能够在较高效率的情况下得出优解；

缺点：需要设定初值，即使是二分法分得再细，也会存在一些不能设置在所有解附近的情况；当然二分法相较随机点法可控，可以采取缩短区间或者提高二分次数的方法缓解这一问题；但是缩短区间意味着你需要大致知道解一定的存在的范围；提高二分次数意味着 2^n+1 二分次数会爆炸😂。

在下抛砖引玉，需要的可以直接使用，有大佬有更合适判断初值的方法，有意望告知，万分感谢！

raymoodlmd

关注

21
点赞
踩
11

收藏

觉得还不错? 一键收藏
1
评论
C# 的多项式方程求解的基本实现

先已知的方法有很多，例如穷举法、二分法、牛顿-莱布利兹法、梯度法、碰撞体积法等等，各有优劣。基本逻辑，当我们设置一个初值进行迭代的时候，最终会获得与该初值距离最近的解；意味着如果方程有多解，我们就必须在每个解的周围都需要至少1个设置的初值。步骤：1、我们要写一个根据已知数据，求拟合方程参数的方法。
复制链接

扫一扫