最小二乘法（高斯）

# include<cstdio>
# include<cmath>
# define N 12   ///12个值
# define M 4
# include<cstring>
# include<iostream>
using namespace std;
///这是数据初始化。
double x[N]= {0,5,10,15,20,25,30,35,40,45,50,55};
double y[N]= {0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64};
double xx[7];
double yy[4];
///这是要用数组，高斯消元要用
double a[5][5]=
{
    {0,0,0,0,0},
    {0,0,0,0,0},
    {0,0,0,0,0},
    {0,0,0,0,0},
    {0,0,0,0,0},
};
double b[5]={0,0,0,0,0};
double nei(double t)
{
    return b[1]*t+b[2]*t*t+b[3]*t*t*t;
}
int main()
{
    int n=3,k,i,j;
    double sum,s,num;

    ///这里是求x[i]、x[i]平方、、、、x[i]六次方的累和。就是正规表达式中的各项要求的x[i]部分。
    for(i=1; i<=6; i++)
    {
        sum=0;
        for(j=0; j<N; j++)
        {
            s=1;
            for(k=1; k<=i; k++)
                s=s*x[j];
            sum=sum+s;
        }
        xx[i]=sum;
    }
    ///这里是看结果对不对。
    for(i=1; i<=6; i++)
    {
        printf("xx[%d]=%lf\n",i,xx[i]);
    }
    ///这里是求正规表达式中的y[i]、、、、x[i]*x[i]*x[i]*y[i]的累和。就是正规表达式中的有y[i]的部分。
    for(i=1; i<=4; i++)
    {
        sum=0;
        for(j=0; j<12; j++)
        {
            s=1;
            for(k=1; k<i; k++)
                s=s*x[j];
            sum=sum+s*y[j];
        }
        yy[i]=sum;
    }
    ///这里是看结果对不对。
    for(i=1; i<=4; i++)
    {
        printf("yy[%d]=%lf\n",i,yy[i]);
    }
    ///这里是把上面求得的x[i],、、、、等放到数组里。
    for(i=1; i<=3; i++)
    {
        b[i]=yy[i];
        for(j=1; j<=3; j++)
        {
            a[i][j]=xx[i+j-1];
        }
    }
     ///这里是看结果在数组放的位置对不对。
    /*
    for(i=1;i<=4;i++)
    {
        for(j=1;j<=4;j++)
            printf("%lf  ", a[i][j]);
        printf("\n");
    }
    */
    /*
    for(i=1;i<=4;i++)
        printf("%lf  ",b[i]);
    printf("\n");
    for(i=1;i<=4;i++)
    {
        for(j=1;j<=4;j++)
            printf("%lf  ", a[i][j]);
        printf("\n");
    }
    */
    ///后面是高斯消元的部分。
    k=1;
    do
    {
        for(i=k+1; i<=n; i++)
            a[i][k]=a[i][k]/a[k][k];

        for(i=k+1; i<=n; i++)
        {
            for(j=k+1; j<=n; j++)
                a[i][j]=a[i][j]-a[i][k]*a[k][j];
            b[i]=b[i]-a[i][k]*b[k];
        }
        if(k==n-1)
            break;
        k++;
    }
    while(1);
    printf("\n");
    for(i=1;i<=4;i++)
    {
        for(j=1;j<=4;j++)
            printf("%lf  ", a[i][j]);
        printf("\n");
    }
    b[n]=b[n]/a[n][n];
    for(i=n-1; i>=1; i--)
    {
        s=0;
        for(j=i+1; j<=n; j++)
            s+=a[i][j]*b[j];
        b[i]=(b[i]-s)/a[i][i];
    }
    ///这里的b[1],b[2],b[3]表示的是拟合函数中的a1,a2,a3.
    for(i=1; i<=n; i++)
        printf("b[%2d]=%lf\n",i,b[i]);
    ///这里输出的是12个值的每一个误差。
    for(i=0;i<12;i++)
    {
        printf("y[%d]=%lf   ",i,y[i]);
        printf("%lf  ",nei(x[i]));
        printf("error  %lf  ",y[i]-nei(x[i]));
        printf("\n");
    }
    printf("\n");
    ///这里求得是12个点平方和的总误差
    num=0;
    for(int i=0;i<12;i++)
    {
        num=(y[i]-nei(x[i]))*(y[i]-nei(x[i]));
    }
    printf("%lf ",num);
    return 0;
}