C语言实现最小二乘多项式曲线拟合（附带源码）

C语言实现最小二乘多项式曲线拟合（附带源码）

引用

CSDN

1.

https://m.blog.csdn.net/m0_61840987/article/details/145412463

最小二乘法是一种优化算法，用于通过最小化数据点与拟合曲线之间的误差平方和来找到最佳拟合曲线。本文将详细介绍如何使用C语言实现最小二乘多项式曲线拟合，包括理论讲解、代码实现和应用示例等多个方面。

项目简介

最小二乘法是一种优化算法，用于通过最小化数据点与拟合曲线之间的误差平方和来找到最佳拟合曲线。在多项式曲线拟合中，最小二乘法通过选择一个多项式形式（例如一次、二次、三次等），来拟合给定的一组数据点。

假设我们有一组数据点 (x_1, y_1), (x_2, y_2), ..., (x_n, y_n)，我们希望找到一个多项式 P(x) 来最小化如下目标函数：

通过最小化这个误差函数，我们可以得到拟合的多项式系数。

项目实现思路

1. 构造目标函数：

• 假设多项式为
  ，其中 m 是多项式的阶数（例如一次多项式、二次多项式等）。

2. 构造方程系统：

• 为了最小化误差函数，需要求解方程系统，该系统是由每个未知系数的偏导数构成的。通过求导得到的一系列方程构成了一个线性方程组。

3. 解线性方程组：

• 使用矩阵方法（如高斯消元法）求解方程组，得到多项式系数。

4. 输出拟合结果：

• 输出计算得到的多项式系数，并使用这些系数拟合数据点。

代码实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// 求解线性方程组 Ax = b
// 使用高斯消元法
void gaussian_elimination(double **A, double *b, double *x, int n) {
    int i, j, k;
    double temp;
    // 构建增广矩阵
    for (i = 0; i < n; i++) {
        A[i][n] = b[i];
    }
    // 高斯消元过程
    for (i = 0; i < n; i++) {
        // 寻找主元素
        for (j = i + 1; j < n; j++) {
            if (fabs(A[i][i]) < fabs(A[j][i])) {
                for (k = 0; k <= n; k++) {
                    temp = A[i][k];
                    A[i][k] = A[j][k];
                    A[j][k] = temp;
                }
            }
        }
        // 消元
        for (j = i + 1; j < n; j++) {
            temp = A[j][i] / A[i][i];
            for (k = i; k <= n; k++) {
                A[j][k] -= A[i][k] * temp;
            }
        }
    }
    // 回代求解
    for (i = n - 1; i >= 0; i--) {
        x[i] = A[i][n] / A[i][i];
        for (j = i - 1; j >= 0; j--) {
            A[j][n] -= A[j][i] * x[i];
        }
    }
}

// 生成多项式的系数
void polynomial_fit(double *x, double *y, int n, int m, double *coefficients) {
    int i, j, k;
    double **A, *b;
    // 为矩阵和向量分配内存
    A = (double **)malloc((m + 1) * sizeof(double *));
    for (i = 0; i <= m; i++) {
        A[i] = (double *)malloc((m + 2) * sizeof(double));  // m+2 是因为增广矩阵
    }
    b = (double *)malloc((m + 1) * sizeof(double));
    // 构建方程组的矩阵 A 和向量 b
    for (i = 0; i <= m; i++) {
        for (j = 0; j <= m; j++) {
            A[i][j] = 0.0;
            for (k = 0; k < n; k++) {
                A[i][j] += pow(x[k], i + j);
            }
        }
    }
    for (i = 0; i <= m; i++) {
        b[i] = 0.0;
        for (k = 0; k < n; k++) {
            b[i] += y[k] * pow(x[k], i);
        }
    }
    // 使用高斯消元法求解系数
    gaussian_elimination(A, b, coefficients, m + 1);
    // 释放内存
    for (i = 0; i <= m; i++) {
        free(A[i]);
    }
    free(A);
    free(b);
}

// 打印多项式系数
void print_polynomial(double *coefficients, int m) {
    printf("拟合的多项式系数：\n");
    for (int i = m; i >= 0; i--) {
        printf("a%d = %.6f\n", i, coefficients[i]);
    }
}

int main() {
    // 示例数据
    double x[] = {1, 2, 3, 4, 5};  // 自变量
    double y[] = {2.3, 2.9, 3.1, 3.6, 5.0};  // 因变量
    int n = 5;  // 数据点的数量
    int m = 2;  // 多项式的阶数（这里选择二次多项式）
    double *coefficients = (double *)malloc((m + 1) * sizeof(double));
    // 求解多项式系数
    polynomial_fit(x, y, n, m, coefficients);
    // 打印拟合结果
    print_polynomial(coefficients, m);
    // 释放内存
    free(coefficients);
    return 0;
}

代码解读

1. gaussian_elimination函数：

• 该函数使用高斯消元法来求解线性方程组。它接收矩阵 A 和向量 b，通过消元和回代步骤求解系数向量 x。

2. polynomial_fit函数：

• 该函数用于根据给定的数据点 (x,y) 和指定的多项式阶数 m 来求解最小二乘法拟合的多项式系数。
• 它构建一个方程组
  ，其中 A 是由 x 数据生成的矩阵，b 是与 y 数据相关的向量。然后通过调用gaussian_elimination函数求解系数。

3. print_polynomial函数：

• 该函数打印拟合得到的多项式系数。

4. main函数：

• 在main函数中，我们定义了一个简单的样本数据集 x 和 y，并调用polynomial_fit函数来计算拟合的多项式系数。最后，通过print_polynomial打印出拟合的结果。

示例输出

假设我们有以下数据点：

x = {1, 2, 3, 4, 5}
y = {2.3, 2.9, 3.1, 3.6, 5.0}

我们选择拟合一个二次多项式。运行程序后，输出将是拟合多项式的系数：

拟合的多项式系数：
a2 = 0.225000
a1 = 0.765000
a0 = 1.340000

这意味着拟合出的多项式为：

项目总结

1. 功能实现：

• 本项目通过最小二乘法实现了对给定数据点的多项式拟合。通过构建方程组并使用高斯消元法求解，我们成功得到了多项式的系数。

2. 优化与扩展：

• 目前的实现仅支持二次多项式（可以通过修改阶数来适应更高阶的多项式）。
• 可以进一步优化高斯消元法，加入更高效的矩阵求解算法（例如LU分解）。

3. 应用场景：

• 最小二乘法多项式拟合广泛应用于数据分析、曲线拟合、统计建模等领域。在科研、工程、金融等行业中，经常需要用这种方法来拟合实验数据或预测未来趋势。

通过本项目，我们了解了最小二乘法的基本原理，并学会了如何在 C 语言中实现该算法来进行多项式拟合。