拟合与逼近

2019独角兽企业重金招聘Python工程师标准>>>

#include "stdafx.h"
#include <stdlib.h>
#include <math.h>
#include "FitApproach.h"

//最佳一致逼近的里米兹方法
void hremz(double a, double b, double p[], int n, double eps)
{
 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;
}

//最小二乘曲线拟合
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;
   }
  }
 }
}

double hremzf(double x)
{
 double y;
 y = exp(x);
 return(y);
}
                                   ----根据《C语言常用算法程序集》整理

 

投票>

转载于:https://my.oschina.net/RapidBird/blog/3522