非线性最小二次拟合C语言

非线性最小二乘拟合问题通常可以用 Levenberg-Marquardt (LM) 算法求解，这是一种基于梯度下降的优化算法。下面是一个使用 C 语言实现的非线性最小二乘拟合示例：

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

#define MAXN 1000  // 最大数据点数
#define MAXM 10    // 最大参数数
#define TOL 1e-10  // 迭代精度

// 待拟合的非线性函数 f(x, p)
double f(double x, double *p, int m) {
    double y = p[0];
    for (int i = 1; i < m; i++) {
        y += p[i] * pow(x, i);
    }
    return y;
}

// 求解雅可比矩阵 J(x, p)
void jacobi(double *J, double x, double *p, int m) {
    for (int i = 0; i < m; i++) {
        J[i] = pow(x, i);
    }
}

// Levenberg-Marquardt 算法
void lm(double *x, double *y, int n, double *p, int m) {
    double J[MAXM][MAXN];
    double JT[MAXN][MAXM];
    double H[MAXM][MAXM];
    double g[MAXM];
    double dp[MAXM];

    double lambda = 0.001;  // 初始阻尼系数
    double prev_err = 1e9;

    for (int iter = 0; iter < 1000; iter++) {
        double err = 0.0;

        // 计算雅可比矩阵 J(x, p) 和误差向量 e(x, p)
        for (int i = 0; i < n; i++) {
            double e = y[i] - f(x[i], p, m);
            jacobi(J[i], x[i], p, m);
            err += e * e;
        }

        if (err < TOL) {
            break;
        }

        // 计算 Hessian 矩阵 H
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                double s = 0;
                for (int k = 0; k < n; k++) {
                    s += J[k][i] * J[k][j];
                }
                H[i][j] = s;
            }
        }

        // 计算梯度向量 g
        for (int i = 0; i < m; i++) {
            double s = 0;
            for (int k = 0; k < n; k++) {
                s += J[k][i] * (y[k] - f(x[k], p, m));
            }
            g[i] = s;
        }

        // Levenberg-Marquardt 步骤
        for (int k = 0; k < m; k++) {
            H[k][k] *= (1 + lambda);
        }

        for (;;) {
            // 解线性方程 H * dp = -g
            for (int i = 0; i < m; i++) {
                for (int j = 0; j < m; j++) {
                    JT[i][j] = H[j][i];
                }
            }
            for (int i = 0; i < m; i++) {
                JT[i][m] = -g[i];
            }
            for (int k = 0; k < m; k++) {
                int pivot = k;
                for (int i = k + 1; i < m; i++) {
                    if (fabs(JT[i][k]) > fabs(JT[pivot][k])) {
                        pivot = i;
                    }
                }
                if (pivot != k) {
                    for (int j = k; j <= m; j++) {
                        double tmp = JT[k][j];
                        JT[k][j] = JT[pivot][j];
                        JT[pivot][j] = tmp;
                    }
                }
                for (int i = k + 1; i < m; i++) {
                    double f = JT[i][k] / JT[k][k];
                    for (int j = k + 1; j <= m; j++) {
                        JT[i][j] -= JT[k][j] * f;
                    }
                    JT[i][k] = 0;
                }
            }
            for (int k = m - 1; k >= 0; k--) {
                double s = JT[k][m];
                for (int j = k + 1; j < m; j++) {
                    s -= JT[k][j] * dp[j];
                }
                dp[k] = s / JT[k][k];
            }

            // 更新参数
            double new_err = 0.0;
            double new_p[MAXM];
            for (int i = 0; i < m; i++) {
                new_p[i] = p[i] + dp[i];
            }
            for (int i = 0; i < n; i++) {
                double e = y[i] - f(x[i], new_p, m);
                new_err += e * e;
            }
            if (new_err < err) {
                lambda /= 10;
                for (int i = 0; i < m; i++) {
                    p[i] = new_p[i];
                }
                prev_err = err;
                break;
            } else {
                lambda *= 10;
                if (lambda > 1e9) {
                    lambda = 1e9;
                    break;
                }
            }
        }

        // 迭代结束条件
        if (prev_err - err < TOL) {
            break;
        }
    }
}

int main() {
    double x[MAXN], y[MAXN];
    int n = 0;

    // 读入数据
    while (scanf("%lf%lf", &x[n], &y[n]) == 2) {
        n++;
    }

    // 拟合二次函数
    double p[MAXM] = {0.0};
    int m = 3;
    lm(x, y, n, p, m);

    // 输出拟合结果
    printf("y = %.3lf + %.3lf*x + %.3lf*x^2\n", p[0], p[1], p[2]);

    return 0;
}

上述代码中，`f(x, p, m)` 函数表示非线性函数的形式，`jacobi(J, x, p, m)` 函数计算雅可比矩阵，`lm(x, y, n, p, m)` 函数使用 Levenberg-Marquardt 算法求解最小二乘拟合参数。在主函数中读入数据并调用 `lm` 函数进行拟合，拟合结果输出为一个二次函数的形式。