最小二乘法.cpp

vC++6.0数值计算方法中“函数拟合最小二乘”

#include<iostream>
#include<math.h>
using namespace std;
void main()
{
	double a[3];
	double x[15]={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8};
	double y[15]={33.4,79.50,122.65,159.05,189.15,241.15,238.65,252.50,267.55,280.50,296.65,301.40,310.40,318.15,325.15};
	void Y(double [],double[],int,int,double[]);
	Y(x,y,15,2,a);
	for(int i=0;i<=2;i++)
	{
		cout<<"a["<<i<<"] = "<<a[i]<<endl;
	}
}
void Y(double x[],double y[],int m ,int n,double a[])
{
	int i,j,t;
	double *c=new double[(n+1)*(n+2)];
	double power(int ,double);
	void fun(double *,int, double[]);
	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]);
	}
	fun(c,n+1,a);
	delete c;
}
void fun(double *c,int n,double x[])
{
	int i,j,k,t;
	double 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;i>=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));
	}
}
double power(int i,double v)
{
	double a=1;
	while(i--)a*=v;
	return a;
}