当前位置:编程学习 > C/C++ >>

[算法分析]微分求积:复化梯形、复化辛浦生

复化梯形
将积分区间[a,b]划分n等分,步长,求积节点,在每个小区间上应用梯形公式
然后将它们累加求和,作为所求积分I的近似值.
记     
式为复化梯形求积公式,下标n表示将区间n等分。
算法流程 
算法代码
[cpp] 
double f(double x){  
  if(x==0) return 1;  
  else return (sin(x)/x);  
}  
  
double FuhuaTixing(int n,double a,double b){  
  double h = (b-a)/n;  
  double x = a;  
  double s = 0;  
  for(int k=0; k< n-1; k++){  
    x += h;  
    s += f(x);  
  }  
  double T = (f(a)+s*2+f(b))*h/2;  
  return T;  
}  
  
int main(){  
  char ans='n';  
  do{  
    cout<<"请输入积分区间(a,b):"<<endl;  
    double a;  
    double b;  
    cin>>a>>b;  
    cout<<"请输入等分份数n: "<<endl;  
    int n; cin>>n;  
    cout<<"由复化梯形公式球的结果:"<<FuhuaTixing(n,a,b)<<endl;  
    cout<<"是否要继续?(y/n)";  
    cin>>ans;  
  }while(ans == 'y');  
  return 0;  
}            
 
复化辛复生
将积分区间[a,b]划分n等分,记子区间的中点为在每个小区间上应用辛普森公式,则有
 
其中
       
记                           
式为复化辛普森求积公式。
 
算法流程
算法代码
[cpp]  
double f(double x){  
  if (x==0)  
     return 1;  
  else return (sin(x)/x);  
}  
  
double Xinfusheng(double a,double b,int n){  
  double h = (b-a)/n;  
  double x = a+1/2*h;  
  double s = 4*f(x);  
  for(int k=1;k<n;k++){  
    x += 1/2*h;  
    s += 4*f(x);  
    x += 1/2*h;  
    s += 2*f(x);  
  }  
  double T=(f(a)+s+f(b))*h/6;  
  return T;  
}  
  
int main(){  
  char ans='n';  
  do{  
    cout<<"请输入积分区间(a,b):"<<endl;  
    double a;  
    double b;  
    cin>>a>>b;  
    cout<<"请输入等分份数n: "<<endl;  
    int n;  
    cin>>n;  
    cout<<"由复化梯形公式球的结果:"<<Xinfusheng(a,b,n)<<endl;  
    cout<<"是否要继续?(y/n)";  
    cin>>ans;  
  }while(ans == 'y');  
  return 0;  
}     
 
实验过程原始记录
分别用复化梯形公式和复化辛浦生公式计算定积分
取n=2,4,8,16,精确解为0.9460831
 
实验结果及分析
1、用复化梯形公式和复化辛甫生公式都能得到较为准确的结果,且等分份数越多,结果的精度越高,梯形公式虽然在等分16份时得到精度与等分4份时相同,但已经越来越接近精确解。辛甫生公式由于C++运算中双精度数值(double)只有7位有效数字的限制,增加等分份数时不容易看出其精度的增加。
2、比较两种方法运算的结果,复化辛甫生公式等分2份时实际要计算5个点的函数值,与复化梯形公式等分4份时计算量基本相同,但得到精度明显复化辛甫生公式要精确很多。
3、复化梯形公式和复化辛甫生公式对于光滑性较差的被积函数都能得到较为精确的结果,而且公式简单,十分利于编译简单的程序由计算机运算,因而得到广泛的应用。
4、实验中的主要误差来自于计算机浮点运算中的截余。
 
补充:软件开发 , C++ ,
CopyRight © 2012 站长网 编程知识问答 www.zzzyk.com All Rights Reserved
部份技术文章来自网络,