切比雪夫（Chebyshev）曲线拟合 （参考《常用算法程序集》） 给定n 个数据点，求切比雪夫意义下的最佳拟合多项式。 函数语句和形参说明： void chir(x,y,n,a,m) double x[n]: 存放给定n 个数据点的X 坐标 double y[n]: 存放给定n 个数据点的Y 坐标 int n: 给定数据点的个数 double a[m+1]: 前 m 个元素返回 m-1 次拟合多项式的m 个系数，最后一个元素返回拟合多项式的偏差最大值。若最后一个元素为负值，则说明在迭代过程中参考偏差不再增大，其绝对值为当前选择的参考偏差。 Int m: 拟合多项式的项数，即拟合多项式的最高次为 m-1。要求m<=n 且 m<=20。若 m>n 或m>20 则自动按 m=min(n-1,20)处理。 源程序： #include "math.h" void chir(x,y,n,a,m) int n,m; double x[],y[],a[]; { 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=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=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; } } } } 例： 取函数f(x)=arctgx 在区间[-1,1]上的101 个点 x(i)=-1.0_0.02i, i=0, 1, … , 100 其相应的函数值为y(i)=f(x(i))。根据此 101 个数据点构造切比雪夫意义下的5 次你和多项式： Ps(x)=a0+a1x+a2x^2+a3x^3+a4x^4+a5x^5 主函数程序如下： #include "math.h" #include "stdio.h" #include "8chir.c" main() { int i; double x[101],y[101],a[7]; for (i=0; i<=100; i++) { x[i]=-1.0+0.02*i; y[i]=atan(x[i]); } chir(x,y,101,a,6); printf("\n"); for (i=0; i<=5; i++) printf("a(%2d)=%e\n",i,a[i]); printf("\n"); printf("MAX(p-f)=%e\n",a[6]); printf("\n"); }
