用C语言写一段高斯牛顿曲线拟合算法

好的，以下是用C语言实现高斯牛顿曲线拟合算法的代码：

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

#define MAXN 1000     // 最大数据点数
#define MAX_ITER 100  // 最大迭代次数
#define EPS 1e-7      // 迭代精度限制

double x[MAXN], y[MAXN]; // 数据点
int n;                   // 数据点数

double func(double a, double b, double c, double x)
{
    return a * exp(-1 * (x - b) * (x - b) / (2 * c * c));
}

// 计算Jacobian矩阵，即一阶导数矩阵
void calculateJacobian(double a, double b, double c, double J[][3])
{
    int i;
    for (i = 0; i < n; i++)
    {
        double t = exp(-1 * (x[i] - b) * (x[i] - b) / (2 * c * c));
        J[i][0] = t;
        J[i][1] = a * t * (x[i] - b) / (c * c);
        J[i][2] = a * t * (x[i] - b) * (x[i] - b) / (c * c * c);
    }
}

// 使用高斯牛顿法拟合
void gaussNewton(double a, double b, double c)
{
    int iter = 0;
    double residuals[MAXN];
    double J[MAXN][3];
    double error = 0.0;

    while (iter < MAX_ITER)
    {
        int i, j;
        double H[3][3] = {0.0};
        double g[3] = {0.0};

        for (i = 0; i < n; i++)
        {
            residuals[i] = y[i] - func(a, b, c, x[i]);
        }

        calculateJacobian(a, b, c, J);

        for (i = 0; i < n; i++)
        {
            for (j = 0; j < 3; j++)
            {
                g[j] += J[i][j] * residuals[i];
            }
        }

        for (i = 0; i < n; i++)
        {
            for (j = 0; j < 3; j++)
            {
                H[j][j] += J[i][j] * J[i][j];
                if (j < 2)
                {
                    H[j][j+1] += J[i][j] * J[i][j+1];
                    H[j+1][j] = H[j][j+1];
                }
            }
        }

        if (fabs(g[0]) < EPS && fabs(g[1]) < EPS && fabs(g[2]) < EPS)
        {
            printf("Converge in %d iteration(s)\n", iter+1);
            break;
        }

        // 求解线性方程组
        double det = H[0][0] * H[1][1] * H[2][2] + H[0][1] * H[1][2] * H[2][0] + H[0][2] * H[1][0] * H[2][1] -
                     H[0][2] * H[1][1] * H[2][0] - H[0][1] * H[1][0] * H[2][2] - H[0][0] * H[1][2] * H[2][1];

        double invH[3][3];
        invH[0][0] = (H[1][1] * H[2][2] - H[1][2] * H[2][1]) / det;
        invH[0][1] = (H[0][2] * H[2][1] - H[0][1] * H[2][2]) / det;
        invH[0][2] = (H[0][1] * H[1][2] - H[0][2] * H[1][1]) / det;
        invH[1][0] = (H[1][2] * H[2][0] - H[1][0] * H[2][2]) / det;
        invH[1][1] = (H[0][0] * H[2][2] - H[0][2] * H[2][0]) / det;
        invH[1][2] = (H[0][2] * H[1][0] - H[0][0] * H[1][2]) / det;
        invH[2][0] = (H[1][0] * H[2][1] - H[1][1] * H[2][0]) / det;
        invH[2][1] = (H[0][1] * H[2][0] - H[0][0] * H[2][1]) / det;
        invH[2][2] = (H[0][0] * H[1][1] - H[0][1] * H[1][0]) / det;

        double delta[3] = {0.0};
        for (i = 0; i < 3; i++)
        {
            for (j = 0; j < 3; j++)
            {
                delta[i] += invH[i][j] * g[j];
            }
        }

        a += delta[0];
        b += delta[1];
        c += delta[2];

        error = 0.0;
        for (i = 0; i < n; i++)
        {
            error += pow(residuals[i], 2);
        }

        printf("%d: a=%lf b=%lf c=%lf error=%lf\n", iter+1, a, b, c, error);

        iter++;
    }
}

int main()
{
    int i;
    printf("请输入数据点数：");
    scanf("%d", &n);
    printf("请输入数据点：\n");
    for (i = 0; i < n; i++)
    {
        scanf("%lf%lf", &x[i], &y[i]);
    }

    double a, b, c;
    printf("请输入初始参数值：");
    scanf("%lf%lf%lf", &a, &b, &c);

    gaussNewton(a, b, c);

    return 0;
}

注意，此代码仅供参考，实际应用中还需要做一些优化和异常处理。