functionapproximation.cpp

能够进行最小平方函数逼近和Chebyshev函数逼近操作。

///////////////////////////////////////////////////
// This program is for bijin.
////////////////////////////////////////////////////
# include <iostream.h>
# include <stdio.h>
# include <math.h>

double LeastsquaresFunctionApproximation(double x[],double y[],int n,double a[],int m,double dt[]);
void ChebyshevFunctionApproximation(double x[],double y[],int n,double a[],int m);

void main()
{
	int i;
	/*
	double x[20],y[20],a[8],dt[3];
	for(i=0;i<20;i++)
	{
		x[i]=0.1*i;
		y[i]=x[i]-exp(-x[i]);
	}
	*/
	double x[5]={0.0,0.25,0.5,0.75,1.0};
	double y[5]={1.0,1.2840,1.6487,2.1170,2.7183};
	double a[5],dt[3];
	double px=0.0;
	for(i=0;i<5;i++)
	{
		a[i]=0.0;
	}
    
	px=LeastsquaresFunctionApproximation(x,y,5,a,5,dt);
	for(i=0;i<5;i++)
	{
		printf("a[%d]=%f \n",i,a[i]);
	}
    printf("x=%f\n",px);
	ChebyshevFunctionApproximation(x,y,5,a,5);
	for(i=0;i<5;i++)
	{
		printf("a[%d]=%f \n",i,a[i]);
	}
    
}


double  LeastsquaresFunctionApproximation(double x[],double y[],int n,double a[],int m,double dt[])
{
    

	int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
	if (m>n)		m=n;
    if (m>20)		m=20;
    for (i=0; i<=m-1; i++)
	    	a[i]=0.0;
    
    z=0.0;
    for (i=0; i<=n-1; i++)
		z=z+x[i]/(1.0*n);
    b[0]=1.0;
	d1=1.0*n;
	p=0.0;
	c=0.0;
    for (i=0; i<=n-1; i++)
	{
		p=p+(x[i]-z);
		c=c+y[i];
	}
    c=c/d1;
	p=p/d1;
    a[0]=c*b[0];

    if (m>1)
    {
		t[1]=1.0;
		t[0]=-p;
        d2=0.0;
		c=0.0;
		g=0.0;
        for (i=0; i<=n-1; i++)
          {
			q=x[i]-z-p; 
			d2=d2+q*q;
            c=c+y[i]*q;
            g=g+(x[i]-z)*q*q;
          }
        c=c/d2;
		p=g/d2;
		q=d2/d1;
        d1=d2;
        a[1]=c*t[1]; 
		a[0]=c*t[0]+a[0];
    }

    for (j=2; j<=m-1; j++)
    {
		s[j]=t[j-1];
		s[j-1]=-p*t[j-1]+t[j-2];
        if (j>=3)
          for (k=j-2; k>=1; k--)
		  { 
              printf("b[%d]=%f\n",k,b[k]);
			  s[k]=-p*t[k]+t[k-1]-q*b[k];
		  }
        s[0]=-p*t[0]-q*b[0];
        d2=0.0;
		c=0.0;
		g=0.0;
        for (i=0; i<=n-1; i++)
        { 
			q=s[j];
            for (k=j-1; k>=0; k--)
              q=q*(x[i]-z)+s[k];
            d2=d2+q*q;
			c=c+y[i]*q;
            g=g+(x[i]-z)*q*q;
        }
        c=c/d2;
		p=g/d2;
		q=d2/d1;
        d1=d2;
        a[j]=c*s[j];
		t[j]=s[j];
        for (k=j-1; k>=0; k--)
        { 
			a[k]=c*s[k]+a[k];
            b[k]=t[k]; 
			t[k]=s[k];
		}
    }

    dt[0]=0.0;
	dt[1]=0.0;
	dt[2]=0.0;
    for (i=0; i<=n-1; i++)
      {
		q=a[m-1];
        for (k=m-2; k>=0; k--)
          q=a[k]+q*(x[i]-z);
        p=q-y[i];
        if (fabs(p)>dt[2]) dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
      }

    double px=0.0;
	for(int ii=0;ii<n;ii++)
	{
		px=px+x[ii];
	}
	px=px/n;
	return(px);
}

void  ChebyshevFunctionApproximation(double x[],double y[],int n,double a[],int m)
 { 
	int m1,i,j,l,ii,k,im,ix[21];
    double h[21],ha,hh,y1,y2,h1,h2,d,hm;
    for (i=0; i<=m; i++) a[i]=0.0;
    if (m>=n) m=n-1;
    if (m>=20) m=19;
    m1=m+1;
    ha=0.0;
    ix[0]=0; ix[m]=n-1;
    l=(n-1)/m; j=l;
    for (i=1; i<=m-1; i++)
      { ix[i]=j; j=j+l;}
    while (1==1)
      { hh=1.0;
        for (i=0; i<=m; i++)
          { a[i]=y[ix[i]]; h[i]=-hh; hh=-hh;}
        for (j=1; j<=m; j++)
          { ii=m1; y2=a[ii-1]; h2=h[ii-1];
            for (i=j; i<=m; i++)
              { d=x[ix[ii-1]]-x[ix[m1-i-1]];
                y1=a[m-i+j-1];
                h1=h[m-i+j-1];
                a[ii-1]=(y2-y1)/d;
                h[ii-1]=(h2-h1)/d;
                ii=m-i+j; y2=y1; h2=h1;
              }
          }
        hh=-a[m]/h[m];
        for (i=0; i<=m; i++)
          a[i]=a[i]+h[i]*hh;
        for (j=1; j<=m-1; j++)
          { ii=m-j; d=x[ix[ii-1]];
            y2=a[ii-1];
            for (k=m1-j; k<=m; k++)
              {
				y1=a[k-1];
			    a[ii-1]=y2-d*y1;
                y2=y1;
				ii=k;
              }
          }
        hm=fabs(hh);
        if (hm<=ha) { a[m]=-hm; return;}
        a[m]=hm; ha=hm; im=ix[0]; h1=hh;
        j=0;
        for (i=0; i<=n-1; i++)
          { if (i==ix[j])
              { if (j<m) j=j+1;}
            else
              { h2=a[m-1];
                for (k=m-2; k>=0; k--)
                  h2=h2*x[i]+a[k];
                h2=h2-y[i];
                if (fabs(h2)>hm)
                  { hm=fabs(h2); h1=h2; im=i;}
              }
          }
        if (im==ix[0]) return;
        i=0;l=1;
        while (l==1)
          { l=0;
            if (im>=ix[i])
              { i=i+1;
                if (i<=m) l=1;
              }
          }
        if (i>m) i=m;
        if (i==(i/2)*2) h2=-hh;
        else h2=hh;
        if (h1*h2>=0.0) ix[i]=im;
        else
          { if (im<ix[0])
              { for (j=m-1; j>=0; j--)
                  ix[j+1]=ix[j];
                ix[0]=im;
              }
            else
              { if (im>ix[m])
                  { for (j=1; j<=m; j++)
                      ix[j-1]=ix[j];
                    ix[m]=im;
                  }
                else ix[i-1]=im;
              }
          }
      }

}