虽然Chebyshev定理从理论上给出了最佳一致逼近的特征性质,但在一般情况下,求最佳一致逼近多项式是困难的,通常只能做近似计算。本文主要介绍Remes算法的利用C++编制的具体程序。
#include <iostream>
#include <cmath>
using namespace std;
void linearfunctions(double x[],double y[],double d[],double u,int n)//解线性方程组的函数
{
int i,j,k;
double **coe=new double* [n+2];
for(i=0; i<n+2; i++)
{
coe[i]=new double[n+2];
}
for(i=0; i<n+2; i++)
{
for(j=0; j<n+1; j++)
{
coe[i][j]=pow(x[i],j);
coe[i][n+1]=pow(-1,i+1);
}
}
for(i=0; i<n+2; i++)
{
if(coe[i][i]==0) break;
}
double *c = new double[n+2];
for(k=0; k<n+1; k++)
{
for(i=k+1; i<n+2; i++)
{
c[i]=coe[i][k]/coe[k][k];
}
for(i=k+1; i<n+2; i++)
{
for(j=0; j<n+2; j++)
{
coe[i][j]=coe[i][j]-c[i]*coe[k][j];
}
y[i]=y[i]-c[i]*y[k];
}
}
if(coe[n+1][n+1]==0) cout<<"Can't solve it!";
double *X=new double [n+2];
X[n+1]=y[n+1]/coe[n+1][n+1];
for(i=n+1; i>=0; i--)
{
double s=0;
for(j=i+1; j<n+2; j++)
{
s=s+coe[i][j]*x[j];
}
X[i]=(y[i]-s)/coe[i][i];
}
for(i=0; i<n+1; i++)
{
d[i]=X[i];
u=fabs(X[n+1]);
}
}
double error(double x[],double y[],double d[],double u,double a,double b,double c,int n)//计算最大模误差的函数
{
double v=0,C,z=0;
for(int i=0; i<n+2; i++)
{
y[i]=sqrt(x[i]); //以具体的f(X)而定,这里是给定的例3.5中的f(X)= √x
}
linearfunctions(x,y,d,u,n);
for(int i=0; i<=10000; i++)
{
C = a + (b-a)*i/10000;
for(int j=0; j<n+1; j++)
{
z+=d[j]*pow(C,j);
}
if(fabs(sqrt(C)-z)>v)
{
v=fabs(sqrt(C)-z);
c=C;
}
}
return v;
}
void change(double x[], double y[], double d[], int n,double a,double b,double c)//调整交错点组的函数
{
double z=0,k=0;
for(int i=0; i<n+1; i++)
{
if(x[i]<c<x[i+1])
{
for(int j=0; j<n+1; j++)
{
z+=d[j]*pow(c,j);
k+=d[j]*pow(x[i],j);
}
if((sqrt(c)-z)*(y[i]-k)>0) x[i]=c;
else x[i+1]=c;
}
}
if(c<x[0])
{
if((sqrt(c)-z)*(y[0]-k)>0) x[0]=c;
else
{
x[0]=c;
for(int i=1; i<n+2; i++)
x[i]=x[i-1];
}
}
if(c>x[n+1])
{
if((sqrt(c-z))*(y[n+1]-k)>0) x[n+1]=(b+x[n+1])/2;
else
{
x[n+1]=c;
for(int i=0; i<n+1; i++)
x[i]=x[i+1];
}
}
}
int main()
{
int i,j,m,n;
double e,a,b,u,v,c;
cout<<"请给定精度:"<<endl;
cin>>e;
cout<<"请给定区间端点值"<<endl;
cin>>a>>b;
cout<<"请给定所需最佳一致逼近多项式的次数:"<<endl;
cin>>n;
double h= (double) (b-a)/(n+1); //n为等分段数,h为步长
double *x= new double[n+2]; //x为交错点组
double *y= new double[n+2]; //y为交错点组对应的函数值
double *d= new double[n+1]; //d为最佳一致逼近多项式的系数
for(i=0; i<n+2; i++)
{
x[i]=a+i*h;
}
v=error(x,y,d,u,a,b,c,n); //v为实际的最大模误差
for(m=0; m<10; m++) //设迭代十次后(即使交错点组仍不符合条件)自动退出
{
cout<<"第"<<m+1<<"次迭代:"<<endl;
change(x,y,d,n,a,b,c);
if(fabs(u-v)<e) // 如果交错点组符合精度要求
{
cout<<"f-pn*近似交错点为:"<<endl;
for(i=0; i<n+2; i++)
{
cout<<x[i]<<" ";
}
cout<<"f(X)的最佳一致逼近多项式的系数从0到n为:"<<endl;
for(i=n; i>=0; i++)
{
cout<<d[i]<<"x^"<<i<<" "; //输出此时的最佳一致逼近多项式
}
break; //跳出循环
}
else if(m<9) //如果交错点组不符合精度要求且迭代次数在十次以内
{
cout<<"此时的交错点组不符合精度要求"<<endl;
}
else
{
cout<<"进行10次迭代后f-pn*近似交错点为:"<<endl;
for(i=0; i<n+2; i++)
{
cout<<x[i]<<" ";
}
cout<<endl<<"f(X)的最佳一致逼近多项式的系数从0到n为:"<<endl;
for(i=n; i>=0; i--)
{
cout<<d[i]<<"x^"<<i<<" "; //输出此时的最佳一致逼近多项式
}
}
}
delete[]x;
delete[]y;
delete[]d;
return 0;
}