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

maleesky

<>积分的源代码,效率很好的
( N/ L" |- Y( ]/ T6 i' {

1 v6 ?0 H) I$ E. _- M6 j//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
: h# P" g4 q* i/ N7 ^' Z5 i//
% }# q! q+ b. `// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. a - Double型变量，积分下限
// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
//
; w4 I& M  D3 i5 b// 返回值：double 型，积分值
  E) A1 h6 E: O& d- G//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
, C3 S8 F" s( j8 m    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
n=1;
8 q, u  g; ]- x% e, r' z" } ep=eps+1.0;
* C  G8 y2 e. ~   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
! Z- s) P5 S4 y* t; _  p=0.0;
; N7 A; D, f' ~3 @6 |; Q        for (i=0;i&lt;=n-1;i++)
% I- }7 B$ M; n        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
3 h4 z. q1 j9 L4 ^6 e6 n+ y3 C        }
$ A' Q9 Y  E" W3 H" F1 {4 G/ R: r2 y        
  p=(y[0]+h*p)/2.0;
! h4 H+ a9 H) z5 A. [        s=1.0;
: }3 U- i$ ]" d+ }( k1 L; `# K        for (k=1;k&lt;=m;k++)
        {
   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
: F: c! N; j! j$ o1 J: |3 K            y[k-1]=p; p=q;
        }</P>
# R/ R+ `; f, z- Q' u+ K<>        ep=fabs(q-y[m-1]);
        m=m+1;
$ s5 [# P3 Q8 f  y[m-1]=q;
$ y% S0 U" x, ?& J, r8 e  n=n+n;
5 _, G3 i4 G: \) N  h=h/2.0;
    }
" Y# L! Y( P% Q   
; `- n( K; f# A* w* [; k return(q);
}</P>