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

maleesky

<>积分的源代码,效率很好的
, l2 a/ i5 g* n. j- @2 {- D, o9 m
5 r. C: M$ t# y
2 Q- @7 M. d6 P; ^. _//////////////////////////////////////////////////////////////////////
& g  n' K- S5 M- b- k// 龙贝格求积法
1 U- ]) c- ~" r//
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
% ~; q* ?2 |4 z1 H// 参数：
8 V' v. t2 v3 m, Q1 V// 1. a - Double型变量，积分下限
// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
( }& \% \9 Y: m" w//
6 U4 i; y1 h+ {// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
% H# `4 K: j2 \. ndouble Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
1 W  R' l' c- h5 w    double y[10],h,ep,p,x,s,q;</P>
, S$ F% H, W8 ^<> // 迭代初值
    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
n=1;
. {% a$ e4 F/ ]' x. a ep=eps+1.0;
   
( t+ ]& k# w; Q2 w // 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
+ a/ T' q" L. H3 u) \    {
' c% i& e0 e7 }/ f8 [- v  p=0.0;
6 Y0 Z$ X  {2 d1 q$ x        for (i=0;i&lt;=n-1;i++)
, k" j7 n$ h3 U7 j( ?, ^7 O( x        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
        }
        
  p=(y[0]+h*p)/2.0;
0 ^/ |7 k" l+ [/ I# Z        s=1.0;
3 r4 U+ r  h5 c8 q! i% R, A, _6 H        for (k=1;k&lt;=m;k++)
        {
   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }</P>
7 U7 s* i- g$ S# c+ g9 A6 F<>        ep=fabs(q-y[m-1]);
4 O: f7 X) N0 ^- g4 W: a  R        m=m+1;
0 o4 z0 `# X  _0 t" _) ]. B  y[m-1]=q;
; g/ ~* S  E: q$ a/ `( j  n=n+n;
  h=h/2.0;
    }
! b- I( z0 C. N/ H% s   
return(q);
( v- |$ N' l1 C9 y}</P>