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

maleesky

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

) m! U' `' Z' x
//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
0 V7 T1 M5 C" C& f9 A1 ?6 T//
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
. G. a5 j; K, t: M4 W// 参数：
) y3 r+ t9 \: }. M// 1. a - Double型变量，积分下限
7 {# G0 C3 M/ a7 k// 2. b - Double型变量，积分上限，要求b&gt;a
// 3. eps - Double型变量，积分精度要求
) v1 ^0 P; W  C' z& b: k//
// 返回值：double 型，积分值
1 t9 j3 W- I  c0 m7 @//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
    int m,n,i,k;
, B' q( [# T. A: W; n    double y[10],h,ep,p,x,s,q;</P>
  ~5 K5 P9 I, v! A" n<> // 迭代初值
1 _. y$ R7 \' R% n; l. Y    h=b-a;
9 K, z/ B3 w. _, ]    y[0]=h*(Func(a)+Func(b))/2.0;
    m=1;
4 }( u3 n% w2 ^, E; Z n=1;
ep=eps+1.0;
   
* P4 T3 w2 I! o3 Y. I // 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
    {
  p=0.0;
        for (i=0;i&lt;=n-1;i++)
        {
   x=a+(i+0.5)*h;
            p=p+Func(x);
5 @8 K/ r- t* ]5 y* {. g- w  I8 b3 J        }
        
  p=(y[0]+h*p)/2.0;
2 Y0 P* O' a2 I9 G' z3 B9 q7 `        s=1.0;
        for (k=1;k&lt;=m;k++)
; X  t7 U. I. \) ?  v# c! u( g- j9 z; B        {
   s=4.0*s;
) Y$ Q+ ?- T# w. E& d  T/ d- U  d            q=(s*p-y[k-1])/(s-1.0);
, I% N+ c5 E: K% {  N2 W            y[k-1]=p; p=q;
        }</P>
$ n1 k6 ^7 r5 B7 P<>        ep=fabs(q-y[m-1]);
        m=m+1;
  y[m-1]=q;
  n=n+n;
  h=h/2.0;
8 [; R2 n" F# @    }
   
5 M/ J( A( J7 V return(q);
% p! H+ j6 [- @/ ^& N1 z}</P>