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

maleesky

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


3 @5 y* |2 B/ J) o4 I: n* ^//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
1 N9 g0 M% O( m  a: ~. V" K; s//
( ?, ^3 ]+ s7 s) D% R8 P// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
) M! V7 R  I, w4 y: F$ T' T$ p// 参数：
# i3 j$ I% z2 p. P// 1. a - Double型变量，积分下限
+ x1 L" X' J" p// 2. b - Double型变量，积分上限，要求b&gt;a
$ S" L& ?% _/ S: y8 n// 3. eps - Double型变量，积分精度要求
//
8 W5 x4 x3 |/ U! N9 Y* m// 返回值：double 型，积分值
  Z( f# z+ K- r; r) [# j//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
5 o2 V& X& w8 o( p. N3 p{
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
0 _, v- r- v7 V    m=1;
/ O5 Z  r% v* ~, Q n=1;
$ x5 U  C* a) [ ep=eps+1.0;
0 \  Y: X. F+ d, D! Z   
: X: m. N/ n9 f$ J7 \5 ~( q // 迭代计算
! J% r3 E7 {; ]: v* O/ b while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
, l' k* p$ u. y: u  p=0.0;
        for (i=0;i&lt;=n-1;i++)
! V5 D! {1 g! }* j% y% l        {
; d0 S1 Y8 u" S7 Q9 u3 b   x=a+(i+0.5)*h;
            p=p+Func(x);
        }
% l0 X8 l" Z( @% F        
5 \6 S. G) ^" T+ X+ P7 X. d  p=(y[0]+h*p)/2.0;
" ^  h2 E2 z  a1 o7 q* S        s=1.0;
        for (k=1;k&lt;=m;k++)
        {
8 q) j# R6 }6 j' o   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
2 K; l5 O4 I* s6 c! h. u2 s* o<>        ep=fabs(q-y[m-1]);
        m=m+1;
4 ?3 g% o" g* }& H0 h2 b7 D' B6 A  y[m-1]=q;
  n=n+n;
0 I6 [: S; }5 d. E  h=h/2.0;
, M: y. [+ x% V- }5 s    }
   
return(q);
}</P>