数学建模社区-数学中国
标题:
龙贝格求积法源代码,效率很好的
[打印本页]
作者:
maleesky
时间:
2005-4-1 16:34
标题:
龙贝格求积法源代码,效率很好的
<
>积分的源代码,效率很好的
# n+ o0 v1 h9 n8 P3 `
' [4 f- Y. r8 H5 r
/ _$ [0 R% L0 u q" v: o+ O
//////////////////////////////////////////////////////////////////////
9 H8 p* Z! L0 V% Q
// 龙贝格求积法
$ d0 z4 v7 T- y A; ~
//
1 h, z; h2 v: `7 j3 ^" u! |- Y
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
& s! G/ E) ?* `4 k4 I. t4 h* q
//
2 V" q0 }/ D4 w$ F4 z2 x. t7 _# b
// 参数：
! F j( U8 Y" A3 x. R$ D0 M
// 1. a - Double型变量，积分下限
% J5 ^# S/ E2 T$ f% t; R) R- n
// 2. b - Double型变量，积分上限，要求b>a
- a3 S, ~( t$ I: d/ u( f
// 3. eps - Double型变量，积分精度要求
1 a& j+ s6 o( X6 C# N& ?+ r
//
' w i, v. R- v5 \7 B# k+ U
// 返回值：double 型，积分值
5 ]* Q2 i* s( l
//////////////////////////////////////////////////////////////////////
8 |! q! p- x L5 V7 y
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
, K' N+ t. Y7 h+ j
{
) E6 X0 P8 G& d5 Y( ^5 @
int m,n,i,k;
$ c% v w* j6 g" v8 J
double y[10],h,ep,p,x,s,q;</P>
" @* |, }! L6 Q+ a+ _. y l4 C
<
> // 迭代初值
8 P7 _; r+ L1 K* ?" J9 T& a6 D
h=b-a;
6 j) p) c( E' c" ]0 w7 [ ~
y[0]=h*(Func(a)+Func(b))/2.0;
; B' [! Z' J) L4 \/ d
m=1;
7 k4 @. K) o7 B
n=1;
, [& k) p- E8 n" R( r7 j9 A5 V
ep=eps+1.0;
& p1 D# |2 i! t
8 K9 f/ V2 y, }7 D: v$ X5 h
// 迭代计算
9 S; R( m) J" }8 N+ y
while ((ep>=eps)&&(m<=9))
! V- ~3 z6 E7 l( N
{
! S7 O# |' C. t i
p=0.0;
/ s! ^7 Z2 x$ O# K
for (i=0;i<=n-1;i++)
0 R4 N. H/ d/ |) }" e/ ]
{
' B3 o _% M1 V$ F1 M# h+ b
x=a+(i+0.5)*h;
7 f0 t5 f2 y+ p
p=p+Func(x);
0 Y1 i8 b7 ]* s m
}
) C( Y3 P6 z/ N# j3 z
3 u# v" R- p; q7 i& A0 d% k
p=(y[0]+h*p)/2.0;
' j! Q$ e0 n* A
s=1.0;
! X# W7 }! ^8 ^ v' R0 j& X* C8 V
for (k=1;k<=m;k++)
7 {. L" \% y9 q1 G. v
{
h5 R6 A# H; M# o, l/ y
s=4.0*s;
% i! x/ O# d$ Z& ^" w
q=(s*p-y[k-1])/(s-1.0);
) B/ W6 `" h" m* O
y[k-1]=p; p=q;
) g8 p( Z' ~) X# |, h
}</P>
1 b( `( d' N' P% W; r& u
<
> ep=fabs(q-y[m-1]);
, l& C8 V5 W; W3 ~3 A6 p1 k! f) S6 c
m=m+1;
3 Q% s! p- A( W2 b5 y
y[m-1]=q;
2 d1 ?1 V, N2 m
n=n+n;
7 C+ O* {' x n' e
h=h/2.0;
, X. C/ g3 }( r( r4 M+ k
}
' D" [! e$ c: `% }9 t9 F- z- d; M
0 T4 `; b+ Q4 X6 `0 q7 D, W
return(q);
1 Z" l x' p4 x' c5 P6 n1 s
}</P>
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
>积分的源代码,效率很好的> // 迭代初值8 P7 _; r+ L1 K* ?" J9 T& a6 D> ep=fabs(q-y[m-1]);