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

maleesky

<>积分的源代码,效率很好的
/ J: V6 s* a# S, v+ l6 D8 J5 a* Z

  U# R: y* S; V/ X7 `% l//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
9 [/ B# X" H5 v) M3 _  ^//
6 e- t1 p: E* @6 s9 x8 B// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
* B, n$ Y; I! @( ]// 参数：
// 1. a - Double型变量，积分下限
! }0 v! S. p. V" U7 s// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
* I" o6 X6 D! Q- v8 V4 b//
// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
) u+ X! v+ m! Odouble Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
& j4 {  Y8 F" ]9 H! i- V( V# g! F2 Q, y    double y[10],h,ep,p,x,s,q;</P>
1 c' r4 b7 j" _. E) x& g/ R7 b<> // 迭代初值
, h' E  c7 c9 L2 s3 D: t; O    h=b-a;
! c4 D! o) p9 |! i& u; A    y[0]=h*(Func(a)+Func(b))/2.0;
3 X6 F9 @+ k" t9 u, K* ?* T% [    m=1;
n=1;
ep=eps+1.0;
   
4 g' l3 P6 Z" F // 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
1 y% }6 S$ _: [. U% t  p=0.0;
        for (i=0;i&lt;=n-1;i++)
, ^" h; v! @4 d5 i- K9 p: U2 I        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
4 f6 e* z' a: O- f        }
( A9 s# s: X8 g2 {        
; t+ E* p" ?# N* c; {  p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
4 ?, z5 o7 |2 D4 v* Z' p! V0 j* |. `        {
" w0 {! \3 d! @' P; R1 R   s=4.0*s;
' P( ^1 H  @0 M% T5 [2 U6 q            q=(s*p-y[k-1])/(s-1.0);
- }# M: {  }* W. O            y[k-1]=p; p=q;
        }</P>
- F7 V% H! q) f- O4 z0 g+ C<>        ep=fabs(q-y[m-1]);
5 m4 i. _2 z+ S4 Y        m=m+1;
$ H' y  K( W2 y9 J: i  y[m-1]=q;
  n=n+n;
, L+ @# i% j1 w8 B1 V1 P& e  h=h/2.0;
    }
  `" |1 b4 q: j, E, s   
return(q);
: l% `3 ]# f8 C* s}</P>