orthrecurs.txt

用递推法产生正交多项式系

/*函数功能：//用递推法产生正交多项式系，即求alpha[j+1]、beta[j];
		入口参数：
			m是数据点数，
			n是拟合的最高阶数，
			float x[],float y[]是对应纵横坐标，
		出口参数：
			a[] 是最小二乘拟合参数
			alpha[]、beta[]是递推系数
注：详细情况请参见《科学计算和C程序集》，蒋长锦编213页
*/

void COkMfcView::orthrecurs(int m, double x[], double y[], double a[], double alpha[], double beta[])
{
	int i,	j,	l ;
	double	d[N],g[N],p,pj,pi,opp,denom,numer,xx ;
	alpha[0] = 0.0 ;
	beta[0] = 0.0 ;
	alpha[1] = 0.0 ;
	p= 0.0 ;

	//计算alpha[1]
	for( i=0 ; i<m ; i++ )
	{
		alpha[1] = alpha[1] + x[i] ;//alpha[1]=1/m{x[1]+x[2]+......+x[n]}
		p = p + 1.0 ;
	}
	alpha[1] = alpha[1] / p ;

	//计算alpha[j+1]、beta[j]
	for( j=1 ; j<N ; j++ )
    {
		denom = 0.0 ;
		numer = 0.0 ;
		for( i=0 ; i<m ; i++ )
	    {
			xx = x[i] ;
			pj = function( j , xx , alpha , beta) ; 
			opp = pj * pj ;
			denom = denom + opp ;
			numer = numer + opp * x[i] ;
		}
		if( denom < ZERO )
	    {
			MessageBox("\n can not determine the orthogonal polynomial fit.	\n please use singular value decomposotion method.\n for minimal least squre solution.\n strike any key to exit.") ;
			goto	hang;
	    }
		beta[j] = denom / p ;
		alpha[j+1] = numer / denom ;
		p = denom ;
	}

	//计算最小二乘拟合参数拟合a[ j ]
	d[0] = 0.0 ;
	g[0] = 0.0 ;
	for( l=0 ; l<m ; l++ )
    {
		d[0] = d[0] + y[l] ;
		g[0] = g[0] + 1.0 ;
    }
	a[0] = d[0] / g[0] ;//cal C0
	for( i=1 ; i<N ; i++ )
	{
		d[i] = 0.0 ;
		g[i] = 0.0 ;
		for( l=0 ; l<m ; l++ )
		{
			xx = x[l] ;
			pi = function(i,xx,alpha,beta) ; 
			g[i] = g[i] + pi * pi ;
			d[i] = d[i] + y[l] * pi ;
	    }
	a[i] = d[i] / g[i] ;
	}

hang:	;
}

//计算各阶正交多项式在X点的函数值
double COkMfcView::function(int ipk, double x, double alpha[], double beta[])
{
	int i ;
	double p[N] ;
	p[0] = 1.0 ;
	p[1] = x - alpha[1] ;
	if( ipk == 1 )
	{
		return(p[1]) ;
	}
	else if( ipk > 1 )
    {
		for( i=1 ; i<ipk ; i++ )
			p[i+1] = ( x - alpha[i+1] ) * p[i] - beta[i] * p[i-1] ;
		return( p[ipk] ) ;
	}
	else
	{

		return 0;
	}
}