ch8.h

/// 科学与工程数值算法 （Java版）

/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH8.h"  拟合与逼近*****/
#ifndef CH8_H_
#define CH8_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"
//*******************************************************************
//作参数的函数指针指向用户自编函数，在具体应用中必须予以确定，函数定义参见文件末
//调用例子: hremz(a,b,p,4,eps,hremzf);直接将函数名作实参即可
//*******************************************************************
void hpir1(double x[],double y[],int n,double a[],int m,double dt[]);//最小二乘曲线拟合
void hchir(double x[],double y[],int n,double a[],int m);//切比雪夫曲线拟合

void hremz(double a,double b,double p[],int n,double eps,
		   double (*hremzf)(double x));//最佳一致逼近的里米兹方法

void hpir2(double x[],double y[],double z[],int n,int m,double a[],int p,int q,double dt[]);//矩形域的最小二乘曲面拟合

//*******************************************************************
void hpir1(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];
    for (i=0; i<=m-1; i++) a[i]=0.0;
    if (m>n) m=n;
    if (m>20) m=20;
    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--)
            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);
      }
    return;
}
/////////////////////////////////////////////////////////////
void hchir(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;
              }
          }
      }
}
/////////////////////////////////////////////////////////////
void hremz(double a,double b,double p[],int n,double eps,
		   double (*hremzf)(double x))
{
    int i,j,k,m;
    double x[21],g[21],d,t,u,s,xx,x0,h,yy;
    if (n>20) n=20;
    m=n+1; d=1.0e+35;
    for (k=0; k<=n; k++)
      { t=cos((n-k)*3.1415926/(1.0*n));
        x[k]=(b+a+(b-a)*t)/2.0;
      }
    while (1==1)
      { u=1.0;
        for (i=0; i<=m-1; i++)
          { p[i]=hremzf(x[i]);
            g[i]=-u; u=-u;
          }
        for (j=0; j<=n-1; j++)
          { k=m; s=p[k-1]; xx=g[k-1];
            for (i=j; i<=n-1; i++)
              { t=p[n-i+j-1]; x0=g[n-i+j-1];
                p[k-1]=(s-t)/(x[k-1]-x[m-i-2]);
                g[k-1]=(xx-x0)/(x[k-1]-x[m-i-2]);
                k=n-i+j; s=t; xx=x0;
              }
          }
        u=-p[m-1]/g[m-1];
        for (i=0; i<=m-1; i++)
          p[i]=p[i]+g[i]*u;
        for (j=1; j<=n-1; j++)
          { k=n-j; h=x[k-1]; s=p[k-1];
            for (i=m-j; i<=n; i++)
              { t=p[i-1]; p[k-1]=s-h*t;
                s=t; k=i;
              }
          }
        p[m-1]=fabs(u); u=p[m-1];
        if (fabs(u-d)<=eps) return;
        d=u; h=0.1*(b-a)/(1.0*n);
        xx=a; x0=a;
        while (x0<=b)
          { s=hremzf(x0); t=p[n-1];
            for (i=n-2; i>=0; i--)
              t=t*x0+p[i];
            s=fabs(s-t);
            if (s>u) { u=s; xx=x0;}
            x0=x0+h;
          }
        s=hremzf(xx); t=p[n-1];
        for (i=n-2; i>=0; i--)
          t=t*xx+p[i];
        yy=s-t; i=1; j=n+1;
        while ((j-i)!=1)
          { k=(i+j)/2;
            if (xx<x[k-1]) j=k;
            else i=k;
          }
        if (xx<x[0])
          { s=hremzf(x[0]); t=p[n-1];
            for (k=n-2; k>=0; k--)
              t=t*x[0]+p[k];
            s=s-t;
            if (s*yy>0.0) x[0]=xx;
            else
              { for (k=n-1; k>=0; k--)
                  x[k+1]=x[k];
                x[0]=xx;
              }
          }
        else
          { if (xx>x[n])
              { s=hremzf(x[n]); t=p[n-1];
                for (k=n-2; k>=0; k--)
                  t=t*x[n]+p[k];
                s=s-t;
                if (s*yy>0.0) x[n]=xx;
                else
                  { for (k=0; k<=n-1; k++)
                      x[k]=x[k+1];
                    x[n]=xx;
                  }
              }
            else
              { i=i-1; j=j-1;
                s=hremzf(x[i]); t=p[n-1];
                for (k=n-2; k>=0; k--)
                  t=t*x[i]+p[k];
                s=s-t;
                if (s*yy>0.0) x[i]=xx;
                else x[j]=xx;
              }
          }
      }
}
/////////////////////////////////////////////////////////////
void hpir2(double x[],double y[],double z[],int n,int m,double a[],int p,int q,double dt[])
{
	int i,j,k,l,kk;
    double apx[20],apy[20],bx[20],by[20],u[20][20];
    double t[20],t1[20],t2[20],xx,yy,d1,d2,g,g1,g2;
    double x2,dd,y1,x1,*v;
    v=(double*)malloc(20*m*sizeof(double));
    for (i=0; i<=p-1; i++)
      { l=i*q;
        for (j=0; j<=q-1; j++) a[l+j]=0.0;
      }
    if (p>n) p=n;
    if (p>20) p=20;
    if (q>m) q=m;
    if (q>20) q=20;
    xx=0.0;
    for (i=0; i<=n-1; i++)
      xx=xx+x[i]/(1.0*n);
    yy=0.0;
    for (i=0; i<=m-1; i++)
      yy=yy+y[i]/(1.0*m);
    d1=1.0*n; apx[0]=0.0;
    for (i=0; i<=n-1; i++)
      apx[0]=apx[0]+x[i]-xx;
    apx[0]=apx[0]/d1;
    for (j=0; j<=m-1; j++)
      { v[j]=0.0;
        for (i=0; i<=n-1; i++)
          v[j]=v[j]+z[i*m+j];
        v[j]=v[j]/d1;
      }
    if (p>1)
      { d2=0.0; apx[1]=0.0;
        for (i=0; i<=n-1; i++)
          { g=x[i]-xx-apx[0];
            d2=d2+g*g;
            apx[1]=apx[1]+(x[i]-xx)*g*g;
          }
        apx[1]=apx[1]/d2;
        bx[1]=d2/d1;
        for (j=0; j<=m-1; j++)
          { v[m+j]=0.0;
            for (i=0; i<=n-1; i++)
              { g=x[i]-xx-apx[0];
                v[m+j]=v[m+j]+z[i*m+j]*g;
              }
            v[m+j]=v[m+j]/d2;
          }
        d1=d2;
      }
    for (k=2; k<=p-1; k++)
      { d2=0.0; apx[k]=0.0;
        for (j=0; j<=m-1; j++) v[k*m+j]=0.0;
        for (i=0; i<=n-1; i++)
          { g1=1.0; g2=x[i]-xx-apx[0];
            for (j=2; j<=k; j++)
              { g=(x[i]-xx-apx[j-1])*g2-bx[j-1]*g1;
                g1=g2; g2=g;
              }
            d2=d2+g*g;
            apx[k]=apx[k]+(x[i]-xx)*g*g;
            for (j=0; j<=m-1; j++)
              v[k*m+j]=v[k*m+j]+z[i*m+j]*g;
          }
        for (j=0; j<=m-1; j++)
          v[k*m+j]=v[k*m+j]/d2;
        apx[k]=apx[k]/d2;
        bx[k]=d2/d1;
        d1=d2;
      }
    d1=m; apy[0]=0.0;
    for (i=0; i<=m-1; i++)
      apy[0]=apy[0]+y[i]-yy;
    apy[0]=apy[0]/d1;
    for (j=0; j<=p-1; j++)
      { u[j][0]=0.0;
        for (i=0; i<=m-1; i++)
	  u[j][0]=u[j][0]+v[j*m+i];
	u[j][0]=u[j][0]/d1;
      }
    if (q>1)
      { d2=0.0; apy[1]=0.0;
        for (i=0; i<=m-1; i++)
          { g=y[i]-yy-apy[0];
            d2=d2+g*g;
            apy[1]=apy[1]+(y[i]-yy)*g*g;
          }
        apy[1]=apy[1]/d2;
        by[1]=d2/d1;
        for (j=0; j<=p-1; j++)
	  { u[j][1]=0.0;
            for (i=0; i<=m-1; i++)
              { g=y[i]-yy-apy[0];
		u[j][1]=u[j][1]+v[j*m+i]*g;
              }
	    u[j][1]=u[j][1]/d2;
          }
        d1=d2;
      }
    for (k=2; k<=q-1; k++)
      { d2=0.0; apy[k]=0.0;
	for (j=0; j<=p-1; j++) u[j][k]=0.0;
        for (i=0; i<=m-1; i++)
          { g1=1.0;
            g2=y[i]-yy-apy[0];
            for (j=2; j<=k; j++)
              { g=(y[i]-yy-apy[j-1])*g2-by[j-1]*g1;
                g1=g2; g2=g;
              }
            d2=d2+g*g;
            apy[k]=apy[k]+(y[i]-yy)*g*g;
            for (j=0; j<=p-1; j++)
	      u[j][k]=u[j][k]+v[j*m+i]*g;
          }
        for (j=0; j<=p-1; j++)
	  u[j][k]=u[j][k]/d2;
        apy[k]=apy[k]/d2;
        by[k]=d2/d1;
        d1=d2;
      }
    v[0]=1.0; v[m]=-apy[0]; v[m+1]=1.0;
    for (i=0; i<=p-1; i++)
      for (j=0; j<=q-1; j++)
        a[i*q+j]=0.0;
    for (i=2; i<=q-1; i++)
      { v[i*m+i]=v[(i-1)*m+(i-1)];
        v[i*m+i-1]=-apy[i-1]*v[(i-1)*m+i-1]+v[(i-1)*m+i-2];
        if (i>=3)
          for (k=i-2; k>=1; k--)
            v[i*m+k]=-apy[i-1]*v[(i-1)*m+k]+
                     v[(i-1)*m+k-1]-by[i-1]*v[(i-2)*m+k];
        v[i*m]=-apy[i-1]*v[(i-1)*m]-by[i-1]*v[(i-2)*m];
      }
    for (i=0; i<=p-1; i++)
      { if (i==0) { t[0]=1.0; t1[0]=1.0;}
        else
          { if (i==1)
              { t[0]=-apx[0]; t[1]=1.0;
                t2[0]=t[0]; t2[1]=t[1];
              }
            else
              { t[i]=t2[i-1];
                t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];
                if (i>=3)
                  for (k=i-2; k>=1; k--)
                    t[k]=-apx[i-1]*t2[k]+t2[k-1]
                         -bx[i-1]*t1[k];
                t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];
                t2[i]=t[i];
                for (k=i-1; k>=0; k--)
                  { t1[k]=t2[k]; t2[k]=t[k];}
              }
          }
        for (j=0; j<=q-1; j++)
          for (k=i; k>=0; k--)
            for (l=j; l>=0; l--)
	      a[k*q+l]=a[k*q+l]+u[i][j]*t[k]*v[j*m+l];
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
    for (i=0; i<=n-1; i++)
      { x1=x[i]-xx;
        for (j=0; j<=m-1; j++)
          { y1=y[j]-yy;
            x2=1.0; dd=0.0;
            for (k=0; k<=p-1; k++)
              { g=a[k*q+q-1];
                for (kk=q-2; kk>=0; kk--)
                  g=g*y1+a[k*q+kk];
                g=g*x2; dd=dd+g; x2=x2*x1;
              }
            dd=dd-z[i*m+j];
            if (fabs(dd)>dt[2]) dt[2]=fabs(dd);
            dt[0]=dt[0]+dd*dd;
            dt[1]=dt[1]+fabs(dd);
          }
      }
    free(v);
    return;
}
#endif
/////////////////////////////////////////////////////////////
/*
double hremzf(double x);
{ 
	double y;
    y=exp(x);
    return(y);
}
/**/