数据拟合终结版.c

计算方法数据拟合

#include<stdio.h>
#include<math.h>
#define MAX_N 20             //定义（x_i,y_i）的最大维数
typedef struct tagPOINT      //定义点结构
{double x;
 double y;
}POINT;
int main()
{
 int n;
 int i,j,k,t,i_k,q;                      //q也是循环变量
 double s,a,b;
 POINT points[MAX_N+1];
 double l[MAX_N+1][MAX_N+1];             //设顶定初始值为零
 double m[MAX_N+1];                      //存储正规方程组的右端项
 double p[2*MAX_N+1]={};                    //用p[k]存储x_i的k次幂的和
 printf("请输入被拟合点的数目,n∈(1,%d)：\n",MAX_N);
 scanf("%d",&n);
 printf("请输入被拟合点（x_i,y_i），i=0,1,2,...,%d:\n",n);
 for(i=0;i<n;i++)
  scanf("%lf%lf",&points[i].x,&points[i].y);
 printf("请输入拟合多项式的次数k(k∈[1,n]):\n");
 scanf("%d",&k);
 for(j=0;j<=2*k;j++)                  //计算x_i的k次幂的和
 {
  p[j]=0;
  for(i=0;i<n;i++)
  {
   p[j]=p[j]+pow(points[i].x,j);
  }
//  printf("%lf\n",p[j]);
 }
 for(j=0;j<=k;j++)                    //计算正规方程组的右端项
 {
  m[j]=0;
  for(i=0;i<n;i++)
  {
   m[j]=m[j]+points[i].y*pow(points[i].x,j);
  }
 }
 for(i=0;i<=k;i++)                    //计算正规方程组的系数矩阵
 {
  for(j=0;j<=k;j++)
  {
   t=i+j;
   l[i][j]=p[t];
  }
 }
 for(i=0;i<=k;i++)                   //输出正规矩阵方程组
 {
  printf("\n");
  for(j=0;j<=k;j++)
  printf("%lf ",l[i][j]);
  printf(" %lf ",m[i]);
 }

//用高斯列主元消去法解正规方程组 ，即求拟合多项式的系数

 for(q=0;q<=k;q++)
 {
  s=0;                      //开始按列选主元
  i_k=q;
  for(i=q;i<=k;i++)
  {
   if(fabs(l[i][q])>s)
   {
    s=fabs(l[i][q]);
    i_k=i;                  //记录主元所在行的指标
   } //if
  }  //for i
  if(s==0)
  {
   printf("系数矩阵奇异！\ndetA=0\n");
   break;
  } //if
  else if(i_k!=q)
  {
   for(j=q;j<=k;j++)       //换行     书上的算法这里错了！！！
   {
    a=l[i_k][j];
    l[i_k][j]=l[q][j];
    l[q][j]=a;
   } //for j
   b=m[i_k];
   m[i_k]=m[q];
   m[q]=b;
  } //else if
  for(i=q+1;i<=k;i++)           //计算乘子
  {
   l[i][q]=l[i][q]/l[q][q];
  } //for i
  for(i=q+1;i<=k;i++)           //消元
  {
   for(j=q+1;j<=k;j++)
   {
    l[i][j]=l[i][j]-l[i][q]*l[q][j];
   } //for j
   m[i]=m[i]-l[i][q]*m[q];
  }  //for i
 } //for q
  for(i=k;i>=0;i--)          //回代；临时用b来存储这里的和
  {
   b=0;
   for(j=i+1;j<=k;j++)
   {b+=l[i][j]*m[j];}
   m[i]=(m[i]-b)/l[i][i];
  } //for i
 printf("\n解为x_i=\n");     //输出正规方程组的解，即拟合多项式的系数
 for(i=0;i<=k;i++)
 printf("%lf\n",m[i]);
 return 0;
}