五阶最小二乘+迭代方法曲线拟合

在最小二乘的基础上，对误差进行拟合，将误差拟合的系数回代，逼近真实值，这个方法是错的，只是按客户要求做一下。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define M 1000
#define N 6
double A[N]= {1.0,2.0,3.0,4.0,5.0,6.0},A_fit[N],A_e_fit[N];
double X[M],Y[M],YY[M],Y_e[M];//原始数据
const char PATH[64]="C:/Users/ren/Desktop/j2453/error/";
/**
 * 说明：5阶拟合函数
 * x：x坐标值
 * y：y坐标值
 * num：数据个数
 * a,b,c,d,e,f：返回拟合系数
 */
void Fit5(double *x, double *y, int num,
          double *a, double *b, double *c, double *d, double *e, double *f)
{
    double sum_x = 0;
    double sum_x2 = 0;
    double sum_x3 = 0;
    double sum_x4 = 0;
    double sum_x5 = 0;
    double sum_x6 = 0;
    double sum_x7 = 0;
    double sum_x8 = 0;
    double sum_x9 = 0;
    double sum_x10 = 0;
    double sum_y = 0;
    double sum_xy = 0;
    double sum_x2y = 0;
    double sum_x3y = 0;
    double sum_x4y = 0;
    double sum_x5y = 0;

    for (int i = 0; i < num; ++i)
    {
        sum_x += x[i];
        sum_x2 += x[i] * x[i];
        sum_x3 += x[i] * x[i] * x[i];
        sum_x4 += x[i] * x[i] * x[i] * x[i];
        sum_x5 += x[i] * x[i] * x[i] * x[i] * x[i];
        sum_x6 += x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
        sum_x7 += x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
        sum_x8 += x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
        sum_x9 += x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
        sum_x10 += x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
        sum_y += y[i];
        sum_xy += x[i] * y[i];
        sum_x2y += x[i] * x[i] * y[i];
        sum_x3y += x[i] * x[i] * x[i] * y[i];
        sum_x4y += x[i] * x[i] * x[i] * x[i] * y[i];
        sum_x5y += x[i] * x[i] * x[i] * x[i] * x[i] * y[i];
    }

    sum_x /= num;
    sum_x2 /= num;
    sum_x3 /= num;
    sum_x4 /= num;
    sum_x5 /= num;
    sum_x6 /= num;
    sum_x7 /= num;
    sum_x8 /= num;
    sum_x9 /= num;
    sum_x10 /= num;
    sum_y /= num;
    sum_xy /= num;
    sum_x2y /= num;
    sum_x3y /= num;
    sum_x4y /= num;
    sum_x5y /= num;

    double m01 = (sum_x8 * sum_x10 - sum_x9 * sum_x9) / (sum_x9 * sum_x10);
    double n01 = (sum_x7 * sum_x10 - sum_x8 * sum_x9) / (sum_x9 * sum_x10);
    double i01 = (sum_x6 * sum_x10 - sum_x7 * sum_x9) / (sum_x9 * sum_x10);
    double j01 = (sum_x5 * sum_x10 - sum_x6 * sum_x9) / (sum_x9 * sum_x10);
    double k01 = (sum_x4 * sum_x10 - sum_x5 * sum_x9) / (sum_x9 * sum_x10);
    double l01 = (sum_x4y * sum_x10 - sum_x5y * sum_x9) / (sum_x9 * sum_x10);

    double m02 = (sum_x7 * sum_x10 - sum_x9 * sum_x8) / (sum_x8 * sum_x10);
    double n02 = (sum_x6 * sum_x10 - sum_x8 * sum_x8) / (sum_x8 * sum_x10);
    double i02 = (sum_x5 * sum_x10 - sum_x7 * sum_x8) / (sum_x8 * sum_x10);
    double j02 = (sum_x4 * sum_x10 - sum_x6 * sum_x8) / (sum_x8 * sum_x10);
    double k02 = (sum_x3 * sum_x10 - sum_x5 * sum_x8) / (sum_x8 * sum_x10);
    double l02 = (sum_x3y * sum_x10 - sum_x5y * sum_x8) / (sum_x8 * sum_x10);

    double m03 = (sum_x6 * sum_x10 - sum_x9 * sum_x7) / (sum_x7 * sum_x10);
    double n03 = (sum_x5 * sum_x10 - sum_x8 * sum_x7) / (sum_x7 * sum_x10);
    double i03 = (sum_x4 * sum_x10 - sum_x7 * sum_x7) / (sum_x7 * sum_x10);
    double j03 = (sum_x3 * sum_x10 - sum_x6 * sum_x7) / (sum_x7 * sum_x10);
    double k03 = (sum_x2 * sum_x10 - sum_x5 * sum_x7) / (sum_x7 * sum_x10);
    double l03 = (sum_x2y * sum_x10 - sum_x5y * sum_x7) / (sum_x7 * sum_x10);

    double m04 = (sum_x5 * sum_x10 - sum_x9 * sum_x6) / (sum_x6 * sum_x10);
    double n04 = (sum_x4 * sum_x10 - sum_x8 * sum_x6) / (sum_x6 * sum_x10);
    double i04 = (sum_x3 * sum_x10 - sum_x7 * sum_x6) / (sum_x6 * sum_x10);
    double j04 = (sum_x2 * sum_x10 - sum_x6 * sum_x6) / (sum_x6 * sum_x10);
    double k04 = (sum_x * sum_x10 - sum_x5 * sum_x6) / (sum_x6 * sum_x10);
    double l04 = (sum_xy * sum_x10 - sum_x5y * sum_x6) / (sum_x6 * sum_x10);

    double m05 = (sum_x4 * sum_x10 - sum_x9 * sum_x5) / (sum_x5 * sum_x10);
    double n05 = (sum_x3 * sum_x10 - sum_x8 * sum_x5) / (sum_x5 * sum_x10);
    double i05 = (sum_x2 * sum_x10 - sum_x7 * sum_x5) / (sum_x5 * sum_x10);
    double j05 = (sum_x * sum_x10 - sum_x6 * sum_x5) / (sum_x5 * sum_x10);
    double k05 = (sum_x10 - sum_x5 * sum_x5) / (sum_x5 * sum_x10);
    double l05 = (sum_y * sum_x10 - sum_x5y * sum_x5) / (sum_x5 * sum_x10);

    double n11 = (m01 * n02 - n01 * m02) / (m01 * m02);
    double i11 = (m01 * i02 - i01 * m02) / (m01 * m02);
    double j11 = (m01 * j02 - j01 * m02) / (m01 * m02);
    double k11 = (m01 * k02 - k01 * m02) / (m01 * m02);
    double l11 = (m01 * l02 - l01 * m02) / (m01 * m02);

    double n12 = (m01 * n03 - n01 * m03) / (m01 * m03);
    double i12 = (m01 * i03 - i01 * m03) / (m01 * m03);
    double j12 = (m01 * j03 - j01 * m03) / (m01 * m03);
    double k12 = (m01 * k03 - k01 * m03) / (m01 * m03);
    double l12 = (m01 * l03 - l01 * m03) / (m01 * m03);

    double n13 = (m01 * n04 - n01 * m04) / (m01 * m04);
    double i13 = (m01 * i04 - i01 * m04) / (m01 * m04);
    double j13 = (m01 * j04 - j01 * m04) / (m01 * m04);
    double k13 = (m01 * k04 - k01 * m04) / (m01 * m04);
    double l13 = (m01 * l04 - l01 * m04) / (m01 * m04);

    double n14 = (m01 * n05 - n01 * m05) / (m01 * m05);
    double i14 = (m01 * i05 - i01 * m05) / (m01 * m05);
    double j14 = (m01 * j05 - j01 * m05) / (m01 * m05);
    double k14 = (m01 * k05 - k01 * m05) / (m01 * m05);
    double l14 = (m01 * l05 - l01 * m05) / (m01 * m05);

    double i21 = (n11 * i12 - i11 * n12) / (n11 * n12);
    double j21 = (n11 * j12 - j11 * n12) / (n11 * n12);
    double k21 = (n11 * k12 - k11 * n12) / (n11 * n12);
    double l21 = (n11 * l12 - l11 * n12) / (n11 * n12);

    double i22 = (n11 * i13 - i11 * n13) / (n11 * n13);
    double j22 = (n11 * j13 - j11 * n13) / (n11 * n13);
    double k22 = (n11 * k13 - k11 * n13) / (n11 * n13);
    double l22 = (n11 * l13 - l11 * n13) / (n11 * n13);

    double i23 = (n11 * i14 - i11 * n14) / (n11 * n14);
    double j23 = (n11 * j14 - j11 * n14) / (n11 * n14);
    double k23 = (n11 * k14 - k11 * n14) / (n11 * n14);
    double l23 = (n11 * l14 - l11 * n14) / (n11 * n14);

    double j31 = (i21 * j22 - j21 * i22) / (i21 * i22);
    double k31 = (i21 * k22 - k21 * i22) / (i21 * i22);
    double l31 = (i21 * l22 - l21 * i22) / (i21 * i22);

    double j32 = (i21 * j23 - j21 * i23) / (i21 * i23);
    double k32 = (i21 * k23 - k21 * i23) / (i21 * i23);
    double l32 = (i21 * l23 - l21 * i23) / (i21 * i23);

    double k33 = (j31 * k32 - k31 * j32) / (j31 * j32);
    double l33 = (j31 * l32 - l31 * j32) / (j31 * j32);

    *f = l33 / k33;
    *e = (l32 - k32 * (*f)) / j32;
    *d = (l23 - k23 * (*f) - j23 * (*e)) / i23;
    *c = (l14 - k14 * (*f) - j14 * (*e) - i14 * (*d)) / n14;
    *b = (l05 - k05 * (*f) - j05 * (*e) - i05 * (*d) - n05 * (*c)) / m05;
    *a = (sum_y - (*f) - sum_x * (*e) - sum_x2 * (*d) - sum_x3 * (*c) - sum_x4 * (*b)) / sum_x5;
}

double Fun_N(double x,const int n,double *a)
{
    double ret=0.0;
    for(int i=0; i<n; i++)
    {
        ret+=(pow(x,i)*a[i]);
    }
    return ret;
}
void save_erro(int x,double err[],int n)
{
    FILE* f;
    char path[64],num[16];
    strcpy(path,PATH);
    itoa(x,num,10);
    strcat(path,num);
    strcat(path,".txt");
    f = fopen(path, "w");
    for(int i=0; i<n; i++)
    {
        fprintf(f,"%f\r",err[i]);
    }
    fclose(f);
}
int main()
{
    //给数据
    for(int i=0; i<M; i++)
    {
        X[i]=0.02*i;
        Y[i]=Fun_N(X[i],N,A);
        YY[i]=Y[i]+0.0003*(double)(rand()%100-50);
    }
    //拟合
    Fit5(X,YY,M,&A_fit[5],&A_fit[4],&A_fit[3],&A_fit[2],&A_fit[1],&A_fit[0]);
    for(int i=0; i<50; i++)
    {
        //求误差相
        for(int j=0; j<M; j++)
        {
            double fit_val=Fun_N(X[j],N,A_fit);
            Y_e[j]=YY[j]-fit_val;

        }
        save_erro(i,Y_e,M);
        //误差行你和
        Fit5(X,Y_e,M,&A_e_fit[5],&A_e_fit[4],&A_e_fit[3],&A_e_fit[2],&A_e_fit[1],&A_e_fit[0]);
        //拟合系数回代
        for(int j=0; j<N; j++)
        {
            A_fit[j]=A_fit[j]+A_e_fit[j];
        }
    }
    return 0;
}