龙贝格求积法源代码,效率很好的

maleesky

<>积分的源代码,效率很好的
; R) K) p0 d" {& ~

//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
; S3 c( H: J7 e' }& u- N8 `) f- z- [//
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
! r3 ^6 e8 ], N. U//
# B0 v! ]  Z9 x! N' Q8 X// 参数：
// 1. a - Double型变量，积分下限
4 w6 Z" m2 i- _// 2. b - Double型变量，积分上限，要求b&gt;a
9 e( I; W) j( V* v4 C% Y: q+ H3 V// 3. eps - Double型变量，积分精度要求
- m- b+ m& h" M' L//
: c! ^1 ?* P% W+ c/ L  Z- u// 返回值：double 型，积分值
" {0 L  P( u! s4 s9 n1 g/ |//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
, D0 i; Y- ?3 ^; c( Z+ V# x4 o: ?& ]    int m,n,i,k;
4 |- }# A, y+ i% i: Y    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
n=1;
ep=eps+1.0;
. y' ^/ d6 T1 B/ L   
0 L! ~$ x6 L+ Z // 迭代计算
9 M. f1 }- v6 G& R7 w. m while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
2 P8 g, @7 Q* x+ ]  y' W    {
2 v8 I/ W+ y( b" P. I' X' A  p=0.0;
2 a5 D( R; ?7 g        for (i=0;i&lt;=n-1;i++)
1 E% }+ X( b3 z        {
  c: @9 `& G. P) g# I6 s' p   x=a+(i+0.5)*h;
            p=p+Func(x);
        }
  |9 ?* F  q) L& `& ], x        
$ a! v8 H/ P* l8 s6 C" W2 h  p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
        {
   s=4.0*s;
) j. }& F: ]5 ?, V  `            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
/ y  Z9 s1 d# P<>        ep=fabs(q-y[m-1]);
        m=m+1;
; ^/ w4 U& T9 P0 z  y[m-1]=q;
  n=n+n;
  h=h/2.0;
    }
   
return(q);
}</P>