最小二乘法（高斯）

# include<cstdio>
# include<cmath>
# define N 12   ///12个值
# define M 4
# include<cstring>
# include<iostream>
using namespace std;
///这是数据初始化。
double x
= {0,5,10,15,20,25,30,35,40,45,50,55};
double y
= {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
=b
/a

;
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;
}