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

maleesky

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


//////////////////////////////////////////////////////////////////////
7 G! |) t' J+ b6 y1 c5 w// 龙贝格求积法
, n  ?! f- f8 ?2 A2 j* s/ H//
4 X3 ^" R1 g2 J// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
& ]2 }) {7 e; O//
// 参数：
# H  `& U5 z5 v6 j: f( S: I// 1. a - Double型变量，积分下限
// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
+ f4 `) d+ m" a0 N7 x0 P. X  i9 k+ y//
// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
7 z, ?. K* m: d  T5 c    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;</P>
  s/ T* m0 A, D8 L<> // 迭代初值
    h=b-a;
( W3 G' e7 Y7 t, f5 B; s    y[0]=h*(Func(a)+Func(b))/2.0;
6 p7 Q  s& |  \2 X+ j% q$ p    m=1;
n=1;
ep=eps+1.0;
   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
  p=0.0;
: J8 Z, Y9 _+ [6 c        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
; K% |4 t  t. ~6 e8 |! v# d- Y" s        }
1 k* l, t( B, q. ~! i        
  p=(y[0]+h*p)/2.0;
        s=1.0;
$ g9 f, t1 \) d8 X  c8 p0 H$ z/ Y! J        for (k=1;k&lt;=m;k++)
2 ]) E& z' E& G0 Y/ Z        {
- j' q8 ~  L+ d/ I5 ~   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
5 E% b9 l. X1 q" r. @<>        ep=fabs(q-y[m-1]);
# \: h0 ^: _6 T% O2 P. x        m=m+1;
" B; f) d: J# B5 l+ M* P, K% s8 \9 [  y[m-1]=q;
1 k8 A6 r5 q: _' L2 M* K  n=n+n;
  h=h/2.0;
' i  Y- o6 y7 e: r" L# m3 E    }
, b! p. M# Y$ H& }3 _2 C   
6 P, ^' J: \- y# F: |. k$ o return(q);
}</P>