虽然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; }
利用Remes算法计算最佳一致逼近多项式
猜你喜欢
转载自blog.csdn.net/weixin_42373618/article/details/80556880
今日推荐
周排行