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

maleesky

<>积分的源代码,效率很好的
& K# w( F5 S$ H. T! q) J7 r

//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
4 ]' ?5 G3 A4 J* w' v3 T* u" P0 E  U//
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
3 K" g6 }! w9 f( t& p; i& \+ Q//
// 参数：
7 [# |) X2 P0 Q( ~/ z$ \// 1. a - Double型变量，积分下限
6 L2 `( @. t! D* T4 e5 e// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
//
// 返回值：double 型，积分值
) u6 s7 o+ R3 J//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
/ R0 ^/ V1 i* z" K4 L    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
- _* N, [4 Q+ I( c    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
6 c. w7 G5 }# J n=1;
ep=eps+1.0;
$ `* \% v& K  I- w4 z" B   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
) O" f+ c! F, p; m/ F' ?; F  p=0.0;
9 G7 S; b4 A$ L% ]        for (i=0;i&lt;=n-1;i++)
& q! \; {  i  j8 Z; D' d$ x' V        {
% h' F8 j" B. o* c   x=a+(i+0.5)*h;
            p=p+Func(x);
1 C0 V) C) ?! U        }
        
  p=(y[0]+h*p)/2.0;
, |& C$ h& ~; h' V) A9 \% o        s=1.0;
        for (k=1;k&lt;=m;k++)
        {
   s=4.0*s;
- D* P" N9 W% ^: r            q=(s*p-y[k-1])/(s-1.0);
# k+ j. [& x$ C! h& J4 j0 W            y[k-1]=p; p=q;
        }</P>
<>        ep=fabs(q-y[m-1]);
        m=m+1;
  y[m-1]=q;
  n=n+n;
  h=h/2.0;
    }
5 Q7 F1 O+ y" G$ M! y  C   
  b4 K4 L4 T, D* r3 T% ] return(q);
1 ~: S  ~2 P: t! k( N6 A2 }$ G. t}</P>