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

maleesky

<>积分的源代码,效率很好的
' X! t7 \; A% D& h
% f2 B; P% u- D7 ?: W
//////////////////////////////////////////////////////////////////////
// 龙贝格求积法
5 ^& ]3 ^5 a3 G3 c//
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
//
. T  {0 a0 G. ~$ D// 参数：
// 1. a - Double型变量，积分下限
" u, r: e- j) Y// 2. b - Double型变量，积分上限，要求b&gt;a
( u7 p4 h0 @! @' U2 d// 3. eps - Double型变量，积分精度要求
2 y7 f( b: V; P% ^# O8 L  X/ Q//
! H/ }$ q9 k0 ^# P4 I8 ~. j// 返回值：double 型，积分值
//////////////////////////////////////////////////////////////////////
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
{
/ d7 ]8 ^0 @" a8 a# K    int m,n,i,k;
6 {( v) t: v! c1 E9 q    double y[10],h,ep,p,x,s,q;</P>
<> // 迭代初值
$ D% t5 B. E! X    h=b-a;
    y[0]=h*(Func(a)+Func(b))/2.0;
: F! a- J. |% w* A9 U( D    m=1;
n=1;
ep=eps+1.0;
   
// 迭代计算
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
6 G7 L$ U) A' [. G1 Z6 ~2 w    {
6 B% x+ r* I9 ~% t& n  p=0.0;
; o* M, t) Q1 U) Q" {        for (i=0;i&lt;=n-1;i++)
/ g# W& p% g8 ~& B" m. Z2 _: q8 N$ `        {
: p% E  e# t1 d   x=a+(i+0.5)*h;
4 c8 H% P" \0 ?            p=p+Func(x);
6 ^9 a: W3 G7 w0 `- F/ B# p) E        }
* r$ J" B9 M; N2 f) _& l3 j& M        
9 _/ r, Y" |% H; T+ @  p=(y[0]+h*p)/2.0;
( e9 q: I( f$ X- W        s=1.0;
( Q6 D1 }" x7 z+ ~        for (k=1;k&lt;=m;k++)
) C# p7 k. g1 [. ~. A        {
% X. z- |# h2 q% j# E9 ~   s=4.0*s;
9 ~9 T* D. N7 E9 y/ `  G            q=(s*p-y[k-1])/(s-1.0);
3 Y' \, G5 a* U$ p" r- {1 u- b            y[k-1]=p; p=q;
        }</P>
<>        ep=fabs(q-y[m-1]);
5 l0 _9 N5 _5 W9 W        m=m+1;
, J8 y2 K3 x* b9 D  y[m-1]=q;
  n=n+n;
' r( o6 d- _& {4 r* [  h=h/2.0;
    }
9 K5 G6 a# D. u9 G   
return(q);
}</P>