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

maleesky

<>积分的源代码,效率很好的
" \/ D# i# n2 R% P7 y% ]. E7 [
( E' ]; U6 R/ q, Q) U: }/ x
//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
( p1 G% U" \- k' T$ ?3 z" @//
) J- Y* C  y, u* x7 X0 Y// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
" G% Q- ~) U! X// 参数：
// 1. a - Double型变量，积分下限
& @" p2 ^3 V! @0 o4 r// 2. b - Double型变量，积分上限，要求b&gt;a
$ V7 G" \8 B% Y% w4 T// 3. eps - Double型变量，积分精度要求
# |8 j: F7 R: ]" ]5 t! G# }# M//
7 G3 x* F3 w" M4 N5 a- \// 返回值：double 型，积分值
- Q3 P- X0 f2 Z) n1 B) q  n7 a//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
- N' K5 h# y* z6 F7 C$ h{
  N7 B! \# Z) U4 ^    int m,n,i,k;
3 V5 K3 R) v. ~& K' q. Y    double y[10],h,ep,p,x,s,q;</P>
) T+ Z" @& ^9 ]% ?$ c. q<> // 迭代初值
    h=b-a;
: J) g# Q5 h( @' Z+ u: f    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
n=1;
ep=eps+1.0;
, i1 p! z: L+ J7 U; D: W3 Z   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
3 S" L, i+ _4 m, W* Q    {
1 p2 p1 x- U" c' Z  ]: R0 v  p=0.0;
: E; o: h8 T- K& h$ \        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
" [1 q+ P* R: O/ l6 A            p=p+Func(x);
        }
        
2 O+ Q9 A# V: e' |  p=(y[0]+h*p)/2.0;
        s=1.0;
2 v  G' X. ^6 o* w6 B; G2 G2 x        for (k=1;k&lt;=m;k++)
; M3 D" p% X7 I  X+ A5 n4 }2 l! B2 p        {
   s=4.0*s;
1 n! M. S! f+ P  |1 i            q=(s*p-y[k-1])/(s-1.0);
" O; G- s. M( e! j) e0 }            y[k-1]=p; p=q;
8 K$ m8 L+ a/ E4 m        }</P>
<>        ep=fabs(q-y[m-1]);
        m=m+1;
& }% A* `, H9 o' d2 @  y[m-1]=q;
  n=n+n;
  h=h/2.0;
3 g  g1 E: Q: ~% z    }
   
return(q);
0 N$ G2 b, m1 T}</P>