c语言实现曲线拟合的最小二乘法

最小二乘法是一种常见的曲线拟合方法，可以用来寻找一条直线或曲线，使其最小化与一组给定数据点之间的误差平方和。以下是C语言实现最小二乘法进行曲线拟合的示例代码。

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

#define MAX_N 100 // 最大数据点数
#define MAX_ORDER 10 // 最大拟合次数

int n, order; // 数据点数、拟合次数
double x[MAX_N], y[MAX_N]; // 数据点
double coef[MAX_ORDER + 1]; // 拟合系数

// 计算幂函数的值
double powi(double x, int n)
{
    double res = 1.0;
    for (int i = 0; i < n; ++i) {
        res *= x;
    }
    return res;
}

// 计算拟合系数
void fit()
{
    double X[MAX_ORDER * 2 + 1] = {0.0}; // 增广矩阵
    double Y[MAX_ORDER + 1] = {0.0}; // 方程右边的向量
    for (int i = 0; i < n; ++i) {
        double xi = x[i], yi = y[i];
        for (int j = 0; j <= order; ++j) {
            X[j] += powi(xi, j) * yi; // 增广矩阵的第一部分
        }
        for (int j = 0; j <= order * 2; ++j) {
            X[j] += powi(xi, j); // 增广矩阵的第二部分
        }
        Y[0] += yi; // 方程右边的向量
        for (int j = 1; j <= order; ++j) {
            Y[j] += powi(xi, j) * yi; // 方程右边的向量
        }
    }
    // 解线性方程组
    for (int k = 0; k <= order; ++k) {
        double akk = X[k];
        if (fabs(akk) < 1e-6) {
            coef[k] = 0.0;
            continue;
        }
        for (int j = k; j <= order; ++j) {
            X[j] /= akk;
        }
        Y[k] /= akk;
        for (int i = k + 1; i <= order; ++i) {
            double aik = X[i * (order * 2 + 1) + k];
            for (int j = k; j <= order * 2; ++j) {
                X[i * (order * 2 + 1) + j] -= aik * X[k * (order * 2 + 1) + j];
            }
            Y[i] -= aik * Y[k];
        }
    }
    for (int k = order; k >= 0; --k) {
        double s = Y[k];
        for (int j = k + 1; j <= order; ++j) {
            s -= X[k * (order * 2 + 1) + j] * coef[j];
        }
        coef[k] = s;
    }
}

// 计算拟合误差
double error()
{
    double err = 0.0;
    for (int i = 0; i < n; ++i) {
        double xi = x[i], yi = y[i];
        double yy = 0.0;
        for (int j = 0; j <= order; ++j) {
            yy += coef[j] * powi(xi, j);
        }
        err += (yy - yi) * (yy - yi);
    }
    return err;
}

int main()
{
    // 读入数据点
    scanf("%d", &n);
    for (int i = 0; i < n; ++i) {
        scanf("%lf%lf", &x[i], &y[i]);
    }
    // 读入拟合次数
    scanf("%d", &order);
    // 计算拟合系数
    fit();
    // 输出结果
    printf("Coefficients:\n");
    for (int i = 0; i <= order; ++i) {
        printf("%.6f ", coef[i]);
    }
    printf("\n");
    printf("Error: %.6f\n", error());
    return 0;
}

该示例代码中，首先读入数据点和拟合次数，然后调用 `fit()` 函数计算拟合系数，最后输出结果和拟合误差。

在 `fit()` 函数中，我们使用增广矩阵法来求解线性方程组，其中增广矩阵的第一部分是 $x_i^j y_i$，第二部分是 $x_i^{j+k}$，方程右边的向量是 $y_i$。然后，在解线性方程组后，我们使用回带法计算拟合系数。

在计算拟合误差时，我们根据公式 $\sum_{i=1}^n (y_i - \hat{y_i})^2$ 来计算误差平方和，其中 $\hat{y_i}$ 是拟合曲线在 $x_i$ 处的值。

注意，该示例代码中使用了一些 C99 标准中的新特性，例如可变长度数组和 `for` 循环中的变量声明。如果您的编译器不支持这些特性，请修改代码以适应您的编译环境。