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

maleesky

<>积分的源代码,效率很好的
' g" n4 ~2 ?" H' T) I! Q4 v  M
1 V6 F: q4 F1 p/ n) n
//////////////////////////////////////////////////////////////////////
: v; n& a* m. Q, x8 P4 t. T// 龙贝格求积法
//
  Q5 K8 M* r1 f+ o$ O, X// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
+ E$ _* J7 z4 t; `& w! [8 t//
// 参数：
, K' G2 f1 f/ k( Z) x// 1. a - Double型变量，积分下限
$ S0 ]# L6 H3 d5 S: C0 ^  e// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
//
// 返回值：double 型，积分值
0 q1 _5 h; b1 z$ I//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
  V+ E7 J( ]2 y  _* V{
0 k$ B6 @$ f: }& p( {$ V# }8 q    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
. x- {# U: b4 l& x, ^    h=b-a;
- g4 S3 Z7 @5 b    y[0]=h*(Func(a)+Func(b))/2.0;
7 @) v' x: h& P5 d4 v$ E    m=1;
n=1;
5 [3 r$ O( t; i3 t ep=eps+1.0;
0 |7 ~& i, G6 W. _/ i. V   
6 w5 j! [8 M$ `9 t- M% a // 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
9 |& u% o2 r* J5 N0 P! J; ?    {
& R/ J1 J) u' U  p=0.0;
( y( F, s) L3 T        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
) I& \# Y# b" Y. e; ?0 c7 b        }
& A& J  n. F6 P: ?( Y0 J        
  p=(y[0]+h*p)/2.0;
0 R4 L3 C/ s" a" k9 i        s=1.0;
        for (k=1;k&lt;=m;k++)
: A# N8 Z; {( o, N  @5 ~        {
/ w. H7 R9 K( L6 ~5 ]: \   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
+ I6 s+ e1 E6 w5 J: G4 `$ T9 `; g<>        ep=fabs(q-y[m-1]);
        m=m+1;
  y[m-1]=q;
: l7 W& Y9 T% f* U8 V* ]0 Q  n=n+n;
  h=h/2.0;
    }
8 t- d4 P- h  E- j   
return(q);
}</P>