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

maleesky

<>积分的源代码,效率很好的
( c# b1 r; E  B7 T7 W
, Y& _- l, v( U1 J
//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
1 @5 p  f7 A2 J$ ^" Q//
1 I- U! f3 i0 W// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
3 ~% u/ k0 r5 a7 G: {7 I//
// 参数：
. W' v6 t" [2 ]7 r* [// 1. a - Double型变量，积分下限
+ \, H- ?5 _% E- {: x9 n- Y" `// 2. b - Double型变量，积分上限，要求b&gt;a
' B5 |2 }6 t  n7 ]. W// 3. eps - Double型变量，积分精度要求
' _& R8 e+ Z1 g, o8 c5 b//
// 返回值：double 型，积分值
0 q7 I& A4 n8 X5 q  z* a; T//////////////////////////////////////////////////////////////////////
, a8 I; t9 m+ f% Wdouble Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
/ [$ b- B6 P; ~0 S, o, }7 S# t* L    double y[10],h,ep,p,x,s,q;</P>
/ K, l! n- g6 y1 P/ x- H( o7 B' V<> // 迭代初值
    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
) N; D  g# {' k    m=1;
n=1;
, X: E- j+ M3 r% H0 z: P ep=eps+1.0;
   
// 迭代计算
( T) O6 g% x  P* ] while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
7 S$ E. i/ A5 j7 m  v  F, {4 k4 }  p=0.0;
        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
, d% p' q( H7 l7 I! _            p=p+Func(x);
        }
: o- Y) i' T! v7 T        
  p=(y[0]+h*p)/2.0;
" d# _3 g  E1 Q6 S6 H        s=1.0;
        for (k=1;k&lt;=m;k++)
5 u" a$ v& _  n& d. ^        {
   s=4.0*s;
3 b* u1 [5 @* e7 n0 x) \            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
1 J7 B6 E8 a' r7 p- W2 g1 v$ C        }</P>
/ S* F2 r' P' r4 O. E) v- n( T<>        ep=fabs(q-y[m-1]);
. Z& N' t' A. x& k1 t) a        m=m+1;
  y[m-1]=q;
6 A' L8 v& t2 x5 c  n=n+n;
& ]- y, B4 V# O# p3 P5 f* j  h=h/2.0;
    }
* K0 Y( K* M+ ?4 V7 y. p   
7 k+ v- Q) [* U( L: X- p7 i return(q);
# H) S# ]# z1 g8 f/ G1 t}</P>