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

maleesky

<>积分的源代码,效率很好的


$ r! @$ p0 i: b/ p' C8 {//////////////////////////////////////////////////////////////////////
8 w5 r4 d4 U* x* h5 c// 龙贝格求积法
9 L% B/ g6 @' q1 l1 \//
5 K) [1 g1 X9 _// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. a - Double型变量，积分下限
4 f  f4 ?. y4 R// 2. b - Double型变量，积分上限，要求b&gt;a
( e0 q+ M0 g1 y) j// 3. eps - Double型变量，积分精度要求
! C0 O  H6 j; s3 r4 X9 O3 @7 C//
// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
$ r2 ]' d/ C1 n' T: ?double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
2 l1 c. W& ~* y( q    int m,n,i,k;
8 T: s) @0 A5 V! M: y    double y[10],h,ep,p,x,s,q;</P>
3 R# V0 ?- v* \/ `5 G* {. V# @) ]<> // 迭代初值
' X3 @: D! K- p7 D# Z! i  \    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
5 f0 ]  ~( n- w/ r" o    m=1;
n=1;
; x; n0 ?3 L! x, }6 x ep=eps+1.0;
) `# v! I, N( C# ^0 m   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
+ Q% m  q2 a! r4 V' {0 ?! l    {
) ^4 {9 e" q1 ?3 R7 O& Y4 v4 Q  p=0.0;
: `0 s% ^# t7 v  }9 L7 ]6 I        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
" P. z5 }$ W" N6 S        }
        
  p=(y[0]+h*p)/2.0;
" x6 l, h* |" ?2 |. p        s=1.0;
        for (k=1;k&lt;=m;k++)
8 }/ b! g% r8 g7 A0 ?- G        {
' C7 F5 i  G6 c* i; g( `4 D$ W- ~   s=4.0*s;
/ _' y, n: ]2 }3 @            q=(s*p-y[k-1])/(s-1.0);
% O6 J) |: k3 S; o            y[k-1]=p; p=q;
3 f$ H( h4 Y* g) f2 [3 b& P        }</P>
<>        ep=fabs(q-y[m-1]);
        m=m+1;
  y[m-1]=q;
' s1 O; U, i4 L  n=n+n;
  h=h/2.0;
    }
   
return(q);
" O0 `) L/ D! _6 }3 @% r}</P>