C#多项式拟合算法实现

以下是一个完整的 C# 多项式曲线拟合实现，采用最小二乘法构建正规方程，并使用列主元高斯消元法求解线性方程组。代码不依赖任何第三方库。

using System;

public static class PolynomialFitter
{
    /// <summary>
    /// 对给定数据点进行多项式拟合，返回从低次到高次排列的系数。
    /// </summary>
    /// <param name="x">自变量数组</param>
    /// <param name="y">因变量数组</param>
    /// <param name="degree">多项式阶数（必须 >= 0）</param>
    /// <returns>系数数组 [a0, a1, a2, ..., an]，对应 a0 + a1*x + a2*x^2 + ...</returns>
    public static double[] Fit(double[] x, double[] y, int degree)
    {
        if (x == null) throw new ArgumentNullException(nameof(x));
        if (y == null) throw new ArgumentNullException(nameof(y));
        if (x.Length != y.Length) throw new ArgumentException("x 和 y 长度必须相同");
        if (degree < 0) throw new ArgumentException("阶数不能为负数");
        if (x.Length < degree + 1) throw new ArgumentException("数据点数必须大于等于阶数+1");

        int n = degree + 1;                     // 未知数个数
        double[,] normalMatrix = new double[n, n];
        double[] normalVector = new double[n];

        // 构建立方程 A^T A a = A^T y
        // 矩阵元素 (i, j) = Σ x_k^(i+j)
        // 右端向量 (i) = Σ y_k * x_k^i
        for (int k = 0; k < x.Length; k++)
        {
            double xi = x[k];
            double yi = y[k];
            double pow = 1.0;

            // 利用递推填充下三角部分，同时累积上三角
            for (int i = 0; i < n; i++)
            {
                normalVector[i] += yi * pow;
                double rowPow = pow;
                for (int j = i; j < n; j++)
                {
                    normalMatrix[i, j] += rowPow;
                    rowPow *= xi;
                }
                pow *= xi;
            }
        }

        // 对称填充下三角
        for (int i = 1; i < n; i++)
            for (int j = 0; j < i; j++)
                normalMatrix[i, j] = normalMatrix[j, i];

        return SolveLinearSystem(normalMatrix, normalVector);
    }

    /// <summary>
    /// 计算多项式在指定 x 处的值。
    /// </summary>
    public static double Predict(double x, double[] coefficients)
    {
        double result = 0.0;
        double pow = 1.0;
        for (int i = 0; i < coefficients.Length; i++)
        {
            result += coefficients[i] * pow;
            pow *= x;
        }
        return result;
    }

    /// <summary>
    /// 使用列主元高斯消元法求解线性方程组 A·x = b
    /// </summary>
    private static double[] SolveLinearSystem(double[,] A, double[] b)
    {
        int n = b.Length;
        double[,] mat = (double[,])A.Clone();
        double[] vec = (double[])b.Clone();

        // 消元过程
        for (int col = 0; col < n; col++)
        {
            // 查找主元
            int maxRow = col;
            double maxVal = Math.Abs(mat[col, col]);
            for (int row = col + 1; row < n; row++)
            {
                double absVal = Math.Abs(mat[row, col]);
                if (absVal > maxVal)
                {
                    maxVal = absVal;
                    maxRow = row;
                }
            }

            if (maxVal < 1e-12)
                throw new InvalidOperationException("正规方程奇异，无法求解。");

            // 交换当前行与最大主元行
            if (maxRow != col)
            {
                for (int j = col; j < n; j++)
                {
                    double tmp = mat[col, j];
                    mat[col, j] = mat[maxRow, j];
                    mat[maxRow, j] = tmp;
                }
                double tmpB = vec[col];
                vec[col] = vec[maxRow];
                vec[maxRow] = tmpB;
            }

            // 消去下方行
            for (int row = col + 1; row < n; row++)
            {
                double factor = mat[row, col] / mat[col, col];
                for (int j = col; j < n; j++)
                    mat[row, j] -= factor * mat[col, j];
                vec[row] -= factor * vec[col];
            }
        }

        // 回代过程
        double[] solution = new double[n];
        for (int row = n - 1; row >= 0; row--)
        {
            double sum = vec[row];
            for (int j = row + 1; j < n; j++)
                sum -= mat[row, j] * solution[j];
            solution[row] = sum / mat[row, row];
        }
        return solution;
    }
}

使用示例

double[] x = { 1, 2, 3, 4, 5 };
double[] y = { 2.1, 3.9, 6.1, 8.2, 9.8 };

// 拟合二次多项式 y = a0 + a1*x + a2*x^2
double[] coeff = PolynomialFitter.Fit(x, y, degree: 2);
Console.WriteLine(string.Join(", ", coeff));  // 输出: a0, a1, a2

// 预测 x = 6 处的值
double predicted = PolynomialFitter.Predict(6, coeff);
Console.WriteLine($"x=6, y≈{predicted:F3}");

算法说明

构建正规方程：对于 m 个数据点 (xₖ, yₖ) 和 n 阶多项式，正规方程矩阵 A（n+1 行 × n+1 列）的元素为 Σ xₖ^(i+j)，右端向量为 Σ yₖ·xₖ^i（i, j 从 0 到 n）。
对称性优化：只计算矩阵上三角，最后复制到下三角以减少计算量。
求解线性方程组：采用列主元高斯消元法保证数值稳定性，当矩阵病态时抛出异常。
系数解释：返回的数组索引对应幂次，即 coeff[0] 是常数项，coeff[1] 是一次项系数，依此类推。

此实现可直接用于 WinForm、WPF 或控制台程序的数值分析场景。

C#多项式拟合算法实现.md

C#多项式拟合算法实现

使用示例

算法说明