基于最小二乘的拟合方法c代码实现空间二次曲面数据拟合

以下是基于最小二乘的拟合方法的 C 代码实现空间二次曲面数据拟合：

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

#define N 10 // 数据点个数
#define M 6  // 参数个数

void least_squares(double x[N], double y[N], double z[N], double A[M][M], double b[M])
{
    int i, j, k;
    double sum_x = 0, sum_y = 0, sum_z = 0, sum_x2 = 0, sum_y2 = 0, sum_xz = 0, sum_yz = 0, sum_x3 = 0, sum_y3 = 0, sum_x2y = 0, sum_xy2 = 0, sum_x4 = 0, sum_y4 = 0, sum_x2z = 0, sum_y2z = 0, sum_xyz = 0;
    for (i = 0; i < N; i++)
    {
        sum_x += x[i];
        sum_y += y[i];
        sum_z += z[i];
        sum_x2 += x[i] * x[i];
        sum_y2 += y[i] * y[i];
        sum_xz += x[i] * z[i];
        sum_yz += y[i] * z[i];
        sum_x3 += x[i] * x[i] * x[i];
        sum_y3 += y[i] * y[i] * y[i];
        sum_x2y += x[i] * x[i] * y[i];
        sum_xy2 += x[i] * y[i] * y[i];
        sum_x4 += x[i] * x[i] * x[i] * x[i];
        sum_y4 += y[i] * y[i] * y[i] * y[i];
        sum_x2z += x[i] * x[i] * z[i];
        sum_y2z += y[i] * y[i] * z[i];
        sum_xyz += x[i] * y[i] * z[i];
    }

    A[0][0] = N;
    A[0][1] = sum_x;
    A[0][2] = sum_y;
    A[0][3] = sum_x2;
    A[0][4] = sum_y2;
    A[0][5] = sum_z;
    A[1][0] = sum_x;
    A[1][1] = sum_x2;
    A[1][2] = sum_x3;
    A[1][3] = sum_x2y;
    A[1][4] = sum_xy2;
    A[1][5] = sum_xz;
    A[2][0] = sum_y;
    A[2][1] = sum_x3;
    A[2][2] = sum_y2;
    A[2][3] = sum_x2y;
    A[2][4] = sum_xy2;
    A[2][5] = sum_yz;
    A[3][0] = sum_x2;
    A[3][1] = sum_x2y;
    A[3][2] = sum_xy2;
    A[3][3] = sum_x4;
    A[3][4] = sum_y2 * sum_x2;
    A[3][5] = sum_x2z;
    A[4][0] = sum_y2;
    A[4][1] = sum_xy2;
    A[4][2] = sum_yz;
    A[4][3] = sum_y2 * sum_x2;
    A[4][4] = sum_x2 * sum_y2;
    A[4][5] = sum_y2z;
    A[5][0] = sum_xz;
    A[5][1] = sum_x2z;
    A[5][2] = sum_y2z;
    A[5][3] = sum_x2z;
    A[5][4] = sum_y2z;
    A[5][5] = sum_x2;

    b[0] = sum_z;
    b[1] = sum_xz;
    b[2] = sum_yz;
    b[3] = sum_x2z;
    b[4] = sum_y2z;
    b[5] = sum_xyz;

    for (i = 0; i < M; i++)
    {
        for (j = i + 1; j < M; j++)
        {
            double temp = A[i][j];
            A[i][j] = A[j][i];
            A[j][i] = temp;
        }
    }

    for (k = 0; k < M; k++)
    {
        double max = A[k][k];
        int max_index = k;
        for (i = k + 1; i < M; i++)
        {
            if (fabs(A[i][k]) > fabs(max))
            {
                max = A[i][k];
                max_index = i;
            }
        }
        if (max_index != k)
        {
            for (j = k; j < M; j++)
            {
                double temp = A[k][j];
                A[k][j] = A[max_index][j];
                A[max_index][j] = temp;
            }
            double temp = b[k];
            b[k] = b[max_index];
            b[max_index] = temp;
        }
        double pivot = A[k][k];
        for (j = k; j < M; j++)
        {
            A[k][j] /= pivot;
        }
        b[k] /= pivot;
        for (i = k + 1; i < M; i++)
        {
            double factor = A[i][k];
            for (j = k + 1; j < M; j++)
            {
                A[i][j] -= factor * A[k][j];
            }
            b[i] -= factor * b[k];
        }
    }

    for (i = M - 1; i >= 0; i--)
    {
        for (j = i - 1; j >= 0; j--)
        {
            b[j] -= A[j][i] * b[i];
            A[j][i] = 0;
        }
    }
}

int main()
{
    double x[N] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    double y[N] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    double z[N] = {1.1, 1.9, 3.2, 4.5, 5.1, 6.3, 7.8, 8.9, 9.8, 10.5};
    double A[M][M];
    double b[M];
    int i, j;

    least_squares(x, y, z, A, b);

    printf("The coefficients of the quadratic surface are:\n");
    printf("z = %f + %fx + %fy + %fx^2 + %fy^2 + %fx*y\n", b[0], b[1], b[2], b[3], b[4], b[5]);

    return 0;
}

这段代码实现了一个简单的空间二次曲面拟合，输入为 $x_i, y_i, z_i$，输出为拟合的二次曲面 $z = a_0 + a_1 x + a_2 y + a_3 x^2 + a_4 y^2 + a_5 x y$ 的系数 $a_0, a_1, a_2, a_3, a_4, a_5$。