最小二乘法.cpp

基于非等距节点的正交多项式的最小二乘法拟合算法之程序实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
const double EPS=1E-10; //运算精度

void readFile(double *&x,double *&y,double *&omiga,int &N)
//读入节点数组x、函数值数组y、权值数组ω及节点数N
{
	FILE *fp;

	fp=fopen(".\\ZXECF.txt","r");
	if (fp==NULL)
	{
		printf("指定位置的文件不存在，请检查！！");
		getch();
		exit(0);
	}

	for (N=1; ;N++)  //计算节点个数N
	{
		fscanf(fp,"%*lf%*lf%*lf");
		if (feof(fp))
			break;
	}

	x=new double[N];
	y=new double[N];
	omiga=new double[N];

	rewind(fp);
	for (int i=0;i<N;i++)
		fscanf(fp,"%lf%lf%lf",&x[i],&y[i],&omiga[i]);

	fclose(fp);
}

void curveFitting(double *x,double *y,double *omiga,double *a,int m,int n)
//x数组存放节点的x[i]值，y数组存放节点的y[i]，omiga数组各节点对应权值
//a数组为拟和的结果多项式的各系数，m为节点个数，n为多项式最高次数
{
	double *alpha,*beta,*Q,*b,*t,*s,*d,*q;
	int i,j,k;
	double sum1,sum2;
	
	alpha=new double[n+1]; //α数组
	beta=new double[n+1];  //β数组
	b=new double[n+1];  //多项式Q(j-1)[x]的系数
	t=new double[n+1];  //多项式Q(j)[x]的系数
	s=new double[n+1];  //多项式Q(j+1)[x]的系数
	Q=new double[n+1];  //临时值d[j]数组
	d=new double[n+1];  //临时值d[j]数组
	q=new double[n+1];  //临时值q[j]数组

	for (i=0;i<=n;i++)  //n次多项式的n+1个系数
		a[i]=0;
	if (n>m) //拟和多项式的最高次数限定不超过m
		n=m;

	//构造Q0(x)
	Q[0]=b[0]=1;  d[0]=sum1=sum2=0;
	for (i=0;i<m;i++)
	{
		sum1+=omiga[i]*x[i]*pow(Q[0],2);
		d[0]+=omiga[i]*pow(Q[0],2);
		sum2+=omiga[i]*y[i]*Q[0];
	}
	q[0]=sum2/d[0];
	alpha[0]=sum1/d[0];
	a[0]=q[0]*b[0];  //计算多项式系数a[0]

	//构造Q1(x)
	t[0]=-alpha[0];
	t[1]=1;
	d[1]=sum1=sum2=0;
	for (i=0;i<m;i++)
	{
		Q[1]=t[0]+t[1]*x[i];
		sum1+=omiga[i]*x[i]*pow(Q[1],2);
		d[1]+=omiga[i]*pow(Q[1],2);
		sum2+=omiga[i]*y[i]*Q[1];
	}
	q[1]=sum2/d[1];
	alpha[1]=sum1/d[1];
	beta[1]=d[1]/d[0];
	a[0]+=q[1]*t[0];
	a[1]=q[1]*t[1];  //计算多项式系数a[1]
	b[1]=t[1];

	//构造Qj(x)
	for (j=2;j<=n;j++)
	{   //计算s[0]--s[j]
		s[j]=t[j-1];
		s[j-1]=-alpha[j-1]*t[j-1]+t[j-2];
		for (k=j-2;k>=1;k--)
			s[k]=-alpha[j-1]*t[k]+t[k-1]-beta[j-1]*b[k];
		s[0]=-alpha[j-1]*t[0]-beta[j-1]*b[0];
		
		sum1=d[j]=sum2=0;
		for (k=0;k<m;k++)
		{
			Q[j]=0;
			for (i=0;i<=j;i++)
				Q[j]+=s[i]*pow(x[k],i);
			sum1+=omiga[k]*x[k]*pow(Q[j],2);
			d[j]+=omiga[k]*pow(Q[j],2);
			sum2+=omiga[k]*y[k]*Q[j];
		}
		q[j]=sum2/d[j];
		alpha[j]=sum1/d[j];
		beta[j]=d[j]/d[j-1];
		for (k=j-1;k>=0;k--)  //多项式拟和
			a[k]+=q[j]*s[k];
		a[j]=q[j]*s[j];
		
		for (k=j-1;k>=0;k--)  //向量赋值：T→B，S→T
		{	b[k]=t[k];
			t[k]=s[k];
		}
		t[j]=s[j];
	}

	//释放动态申请的内存空间
	delete []alpha;
	delete []beta;
	delete []s;
	delete []t;
	delete []b;
	delete []Q;
	delete []d;
	delete []q;
}

void main()
{
	double *x,*y; //x数组存放节点的x[i]值，y数组存放节点的y[i]
	double *omiga; //各节点对应权值数组
	double *a;  //a数组存放n次多项式的n+1个系数
	int m,n,i;  //m为节点个数，编号为0...m-1；n为多项式拟合次数

	readFile(x,y,omiga,m); //从文件中读取数据
	printf("原数据为：\n");
	printf("%-10s%-10s%-10s\n","x","y","节点权值");
	for (i=0;i<m;i++)
		printf("%-10g%-10g%-10g\n",x[i],y[i],omiga[i]);
	printf("\n节点个数为：m=%d\n\n",m);

	n=3;  //3次多项式拟合
	a=new double[n+1]; 

	curveFitting(x,y,omiga,a,m,n);

	printf("%d次多项式拟和结果为：\n",n);
	printf("f(x)=%lg*x^%d",a[n],n);
	for (i=n-1;i>=0;i--)
	{
		if (fabs(a[i])<1E-10)
			continue;
		if (a[i]>0)
			printf("+");
		if (i>1)
			printf("%lg*x^%d",a[i],i);
		else if (i==1)
			printf("%lg*x",a[i]);
		else if (i==0)
			printf("%lg",a[i]);
	}
	printf("\n\n");

	delete []x;
	delete []y;
	delete []omiga;
	delete []a;
}