最小二乘拟合C语言程序

输入文件zuixiaoercheng_input.txt中内容：

5 3

-2 -1 0 1 2

-0.1 0.1 0.4 0.9 1.6

代码部分：

#include

#include

#include

#include

#define p(k) (p+(k)*m)

#define pkx(k) (pkx+(k)*n)

double *x,*f,*p,*a,*b,*d,*pkx,*sx,*xishu;

int m,n;

void zhengjiaonihe();

void size(int *m,int *n);

void input(double x[],double f[],int m,int n);

void output(int n,double pkx[],double xishu[],double R);

double NeiJi(double *F,double *G,double *W);

void GetPk(int k);

void GetXiShu(int k);

double Wucha(int k);

void main()

{

zhengjiaonihe();

system("pause");

}

void zhengjiaonihe()

{

int k,i,j;

double t,R;

printf("\n本程序是用正交函数作最小二乘拟合\n");

size(&m,&n);

x=(double *)malloc(sizeof(double)*m);

f=(double *)malloc(sizeof(double)*m);

p=(double *)malloc(sizeof(double)*m*n);

a=(double *)malloc(sizeof(double)*n);

b=(double *)malloc(sizeof(double)*n);

d=(double *)malloc(sizeof(double)*n);

pkx=(double *)malloc(sizeof(double)*n*n);

sx=(double *)malloc(sizeof(double)*m);

xishu=(double *)malloc(sizeof(double)*n);

input(x,f,m,n);

for (i=0;i

{

for (j=0;j

pkx(i)[j]=0; // pkx 初始化

}

for (k=0;k

{

p(0)[k]=1; // P0(x)=1;

}

pkx(0)[0]=1;

d[0]=NeiJi(f,p(0),p(0))/NeiJi(p(0),p(0),p(0));

GetXiShu(1);

R=Wucha(1);

a[1]=NeiJi(x,p(0),p(0))/NeiJi(p(0),p(0),p(0)); //α1=(xP1,P1)/(P1,P1)

pkx(1)[0]=-a[1];pkx(1)[1]=1;

for (k=0;k

p(1)[k]=x[k]-a[1]; // P1(x)=(x-α1)P0(x) P0(x)=1

d[1]=NeiJi(f,p(1),p(0))/NeiJi(p(1),p(1),p(0));

GetXiShu(2);

R=Wucha(2);

for (k=1;k

{

t=NeiJi(p(k),p(k),p(0));

a[k+1]=NeiJi(p(k),p(k),x)/t; // αk+1=(xPk,Pk)/(Pk,Pk)

b[k]=t/NeiJi(p(k-1),p(k-1),p(0)); //βk=(Pk,Pk)/(Pk-1,Pk-1)

GetPk(k+1); // 计算Pk(x)在各节点的函数值

d[k+1]=NeiJi(f,p(k+1),p(0))/NeiJi(p(k+1),p(k+1),p(0));

// ak=(f,Pk)/(Pk,Pk)

pkx(k+1)[0]=-a[k+1]*pkx(k)[0]-b[k]*pkx(k-1)[0];

for (i=1;i

{ // Pk+1=(x-αk+1)Pk(x)-βkPk-1(x)

pkx(k+1)[i]=pkx(k)[i-1]-a[k+1]*pkx(k)[i]-b[k]*pkx(k-1)[i];

}

GetXiShu(k+1); // 求拟合多项式的系数

R=Wucha(k+1); //计算平方误差

}

GetXiShu(n); // 求拟合多项式的系数

R=Wucha(n); //计算平方误差

output(n,pkx,xishu,R);

}

double NeiJi(double *F,double *G,double *W) //内积函数

{

double s=0;

int i;

for (i=0;i

s+=F[i]*G[i]*W[i];

return s;

} // NeiJi end

void GetPk(int k) // 求正交多项式Pk(x)在节点 x(k) 的值

{

int i;

for (i=0;i

p

(k)[i]=(x[i]-a[k])*p(k-1)[i]-b[k-1]*p(k-2)[i];

}

void GetXiShu(int k) // 求拟合多项式的系数

{

int i,j;

double t;

for (i=0;i

{

t=0;

for (j=0;j

{

t+=d[j]*pkx(j)[i];

}

xishu[i]=t;

}

} // GetXiShu end

double Wucha(int k) // 计算平方误差

{

double r1;

int i,j;

for (i=0;i

{

sx[i]=0;

for (j=0;j

sx[i]=sx[i]+xishu[j]*pow(x[i],j); // 计算拟合函数在各节点的值

}

r1=fabs(NeiJi(f,f,p(0))-NeiJi(f,sx,p(0))); // 计算平方误差

return r1;

} // Wucha end

void input(double x[],double f[],int m,int n)

{ FILE *fp;

char *filename="zuixiaoercheng_input.txt";

int j;

float tmp;

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",&m);

fscanf(fp,"%d",&n);

printf("结点个数m=%d,次数n=%d\n",m,n);

printf("x=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

x[j]=tmp; /* 读取x */

printf("%f ",tmp);

}

printf("\n");

printf("f=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

f[j]=tmp; /* 读取y */

printf("%f ",tmp);

}

fclose(fp);

}

void output(int n,double pkx[],double xishu[],double R)

{ FILE *fp;

char *filename="zuixiaoercheng_output.txt";

int k,i;

if((fp=fopen(filename,"w"))==NULL)

{printf("can not open file.\n");

system("pause");

}

printf("\n正交多项式Pk(x)为: \n");

for (k=0;k

{

{printf("P%d(x)=%lG",k,pkx(k)[0]);

fprintf(fp,"P%d(x)=%lG",k,pkx(k)[0]);

for (i=1;i

if(pkx(k)[i]>=0)

{printf("+%lGx^%d",pkx(k)[i],i);

fprintf(fp,"+%lGx^%d",pkx(k)[i],i);

}

else

{printf("%lGx^%d\n",pkx(k)[i],i);

fprintf(fp,"%lGx^%d\n",pkx(k)[i],i);

}

}

printf("\n");

fprintf(fp,"\n");

}

printf("\n平方误差为: %lf\n",R);

printf("\n拟合多项式为: \n");

printf("F(x)=%lG",xishu[0]);

fprintf(fp,"\n平方误差为: %lf\n",R);

fprintf(fp,"\n拟合多项式为: \n");

fprintf(fp,"F(x)=%lG",xishu[0]);

for (k=1;k

{

if(xishu[k]>=0)

{printf("+%lGx^%d",xishu[k],k);

fprintf(fp,"+%lGx^%d",xishu[k],k);

}

else

{printf("%lGx^%d",xishu[k],k);

fprintf(fp,"%lGx^%d",xishu[k],k);

}

}

printf("\n");

fclose(fp);

}

void size(int *m,int *n)

{FILE *fp;

char *filename="zuixiaoercheng_input.txt";

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",m);

fscanf(fp,"%

d",n);

fclose(fp);

}