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

maleesky

<>积分的源代码,效率很好的

: f3 y/ F) S( i0 _
/ C, k" N) c' u1 h//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
* k' D! Y) M, m0 k* i//
5 x" d4 u/ U; ^7 }  b  I* F// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
7 }! t! Z, @1 C6 j// 参数：
7 G& T/ ?8 _9 L/ e7 a5 r+ A// 1. a - Double型变量，积分下限
2 t% K' v7 p! t0 t8 }. o// 2. b - Double型变量，积分上限，要求b&gt;a
% A4 j' d5 O- P; F' R// 3. eps - Double型变量，积分精度要求
//
// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
% Y$ t) M5 K5 K+ ~, U; Tdouble Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
; m  ~$ ?/ e" ^# M1 ?, V2 W; S    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
- F" M- e% N8 X3 y    m=1;
1 i0 z( Q: F4 ~- s; ?: `5 Z* | n=1;
0 o7 F- n3 T3 }% e  u+ i( ` ep=eps+1.0;
/ i  v( k0 e  g4 S   
// 迭代计算
3 W6 i  v0 w" R* J while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
* D% n9 ]! j* |2 W+ J) R  p=0.0;
8 g/ |' n6 ]/ `. }# F        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
        }
        
  p=(y[0]+h*p)/2.0;
        s=1.0;
& Y" {6 O. m5 X* n, x        for (k=1;k&lt;=m;k++)
        {
   s=4.0*s;
% Q# v* ^2 z+ L" l& z* T            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
) R# J% k* X1 U4 L" I! e' F<>        ep=fabs(q-y[m-1]);
) V! H4 a8 O, a3 h8 R        m=m+1;
2 X/ @) y0 P8 |$ U* h6 d7 Y. B) h  y[m-1]=q;
  n=n+n;
! l( J" c* w7 h% K$ l! k  h=h/2.0;
    }
   
return(q);
3 j$ f2 o3 ]% N}</P>