curvefitting.cpp

基于最小二乘原理的曲线拟合程序-A program of curve fitting based on least squares algorithm.

#include<stdio.h>
#include<math.h>

void main()
{
	int i;
	float a[3];
	float x[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
	float y[16]={4,6.4,8,8.8,9.22,9.5,9.7,9.86,10,10.2,10.32,10.42,10.50,10.55,10.58,10.60};
	void Approx(float[],float[],int,int,float[]);
	Approx(x,y,16,2,a);
	for(i=0;i<=2;i++)
		printf("a[%d]=%f\n",i,a[i]);
}
void Approx(float x[],float y[],int m,int n,float a[])
{
	int i,j,t;
	float *c=new float[(n+1)*(n+2)];
	float power(int,float);
	void ColPivot(float *,int,float[]);
	for(i=0;i<=n;i++)
	{
		for(j=0;j<=n;j++)
		{
			*(c+i*(n+2)+j)=0;
			for(t=0;t<=m-1;t++)
				*(c+i*(n+2)+j)+=power(i+j,x[t]);
		}
		*(c+i*(n+2)+n+1)=0;
		for(j=0;j<=m-1;j++)
			*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);
	}
	ColPivot(c,n+1,a);
	delete c;
}
void ColPivot(float *c,int n,float x[])
{
	int i,j,t,k;
	float p;
	for(i=0;i<=n-2;i++)
	{
		k=i;
		for(j=i+1;j<=n-1;j++)
			if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) k=j;
		if(k!=i)
			for(j=i;j<=n;j++)
			{
				p=*(c+i*(n+1)+j);
				*(c+i*(n+1)+j)=*(c+k*(n+1)+j);
				*(c+k*(n+1)+j)=p;
			}
		for(j=i+1;j<=n-1;j++)
		{
			p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));
			for(t=i;t<=n;t++)
				*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));
		}
	}
	for(i=n-1;j>=0;i--)
	{
		for(j=n-1;j>=i+1;j--)
			(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));
		x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));
	}
}

float power(int i,float v)
{
	float a=1;
	while(i--)a*=v;
	return a;
}