Tauruser

Enjoy Every Day
posts - 34, comments - 95, trackbacks - 0, articles - 5
  C++博客 :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理

Romberg求积(例程)

Posted on 2006-05-29 09:25 Tauruser 阅读(1385) 评论(1)  编辑 收藏 引用 所属分类: 数值计算

先介绍一下Romberg求积。
6.3 外推原理与Romberg求积

6.3.1 复合梯形公式递推化与节点加密

  在计算机上用等距节点求积公式时,若精度不够可以逐步加密节点.设将区间分为n等分,节点,在区间上梯形公式为
             
若节点加密一倍,区间长为,记中点为在同一区间上的复合梯形公式惟
             
于是
       (6.3.1)
它表明是在的基础上再加新节点的函数值之和乘新区间长,而不必用(6.2.6)重新计算,这时有误差估计式
             
,则得
            (6.3.2)
它表明用,其误差近似.这也是在计算机上估计梯形公式误差的近似表达式.若(给定精度),则.
  若在区间[a,b]中做2n等分时,在上用Simpson公式计算,则由(6.2.8)可知
   
它恰好是(6.3.2)中I(f)的近似值,即

它表明用(6.3.2)计算I(f),其精度已由提高到如果再将区间分半,使分为4个小区间,长度为,则可由(6.3.1)计算出,利用复合公式余项(6.2.9)得
         
         
如果,则有
               (6.3.3)
从而有复合Simpson公式的误差估计
         
如果用(6.3.3)近似,即
               (6.3.4)
则精度可达到.类似做法还可继续下去.这样对区间逐次分半,利用公式(6.3.1)逐次递推.再由(6.3.2),(6.3.3)逐次构造出精度愈来愈高的计算积分I(f)的公式,这就是Romberg求积的基本思想.

以下为我自己写的求积程序。

//  RombergIntegral.cpp : 定义控制台应用程序的入口点。
//

#include 
< cmath >
#include 
< iostream >
#include 
< vector >
using   namespace  std;
const   double  PRECISION(. 000001 ); // 精度控制
const  unsigned  int  MAXK( 20 ); // 求解步骤控制
double  RombergIntegral( double  ( * f)( double  x), double  a,  double  b);
vector
< vector < double >>  T; // 用于存储T表
double  f( double  x) // 要求的积分函数
{
    
return  x * sin(x);
}

int  _tmain( int  argc, _TCHAR *  argv[])
{
    cout
<< " 本程序用于求解函数f(x)=x*sin(x)在0到6.28的积分 " << endl;
    cout
<< " 积分结果为: " << RombergIntegral(f, 0 , 6.28 ) << endl;
    cout
<< " 精度为 " << PRECISION << endl;
    
return   0 ;
}


double  RombergIntegral( double  ( * f)( double  x), double  a,  double  b)
{
    
int  k( 0 );
    
double  h = b - a;
    vector
< double >  temp;
    T.push_back(temp);
    T[
0 ].push_back(h * (( * f)(a) + ( * f)(b)) / 2 );
    
for (k = 1 ; 1 ; ++ k)
    
{
        T.push_back(temp);
        T[
0 ].push_back( 0.5 * T[ 0 ][k - 1 ]);
        
for ( int  i = 0 ;i < pow( 2 .,k - 1 ); ++ i)
        
{
            T[
0 ][k] += 0.5 * h * (( * f)(a + h / 2 + i * h));
        }

        
for ( int  i = 1 ;i <= k; ++ i)
            T[i].push_back((pow(
4 .,i) * T[i - 1 ].back() - T[i - 1 ][T[i - 1 ].size() - 2 ]) / (pow( 4 .,i) - 1 ));
        h
/= 2 ;
        
double  temp = T[k].back();
        
if (fabs(T[k].front() - T[k - 1 ].front()) < PRECISION  ||   k == MAXK)  break ; //
    }

    
    
return  T[k].back();
}


// 以上程序在vs2005+win2003下编译运行通过。

Feedback

# re: Romberg求积(例程)  回复  更多评论   

2006-08-06 19:12 by 天涯草
你好厉害啊.佩服!

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理