数学建模社区-数学中国
标题:
龙贝格求积法源代码,效率很好的
[打印本页]
作者:
maleesky
时间:
2005-4-1 16:34
标题:
龙贝格求积法源代码,效率很好的
<
>积分的源代码,效率很好的
3 @" R+ d3 S: w0 T
/ ^) n e B4 Y; ~* ]
, M% w8 F6 w5 ]/ a
//////////////////////////////////////////////////////////////////////
* X2 D4 ~ I7 m9 K" G; S$ s5 \8 c
// 龙贝格求积法
$ m0 b0 E9 k$ R" s, s
//
$ w3 s" D% e: l% A6 J( B+ q' D! m
// 调用时，须覆盖计算函数f(x)值的虚函数double Func(double x)
& E7 I7 i. N- x; M T& s
//
& m6 b% h+ j/ Q' b5 l1 m* \1 ]+ L
// 参数：
* e0 X( M' w& M A# I# r5 ~# _+ H
// 1. a - Double型变量，积分下限
, k: S: h4 ~% J
// 2. b - Double型变量，积分上限，要求b>a
/ a- S I: d! o* Q, S0 V
// 3. eps - Double型变量，积分精度要求
# y7 D% v( z S) H: e, P) R
//
, r8 j9 Z3 H4 x
// 返回值：double 型，积分值
5 S: K% C* W3 v6 \
//////////////////////////////////////////////////////////////////////
, q- E) `" x" D% K7 K: X
double Integral_Romberg(double a, double b, double eps /*= 0.000001*/)
* g8 C$ t# Z. I. N$ C, @
{
# I3 k6 r- M N* ~. j- \! U
int m,n,i,k;
4 O0 i2 Z$ g( w1 S
double y[10],h,ep,p,x,s,q;</P>
5 R, j4 N0 K* b' J M
<
> // 迭代初值
- r y7 Q1 i" D6 r' t
h=b-a;
3 w" @7 a. W6 u d
y[0]=h*(Func(a)+Func(b))/2.0;
: g' e' P. t/ u \, I& \- t" J
m=1;
2 V- C7 v1 P6 i2 ?9 V9 M
n=1;
, }$ S* p$ c9 f- B% \3 z: R
ep=eps+1.0;
( z5 R R4 M+ P/ J* K: |; ^; L' J
7 z4 v$ g0 ]3 {3 U' l
// 迭代计算
; Y( \2 S+ g# Z1 z5 }( E$ U
while ((ep>=eps)&&(m<=9))
4 F P* U; h* k0 N1 i3 d; R
{
) i# p9 @8 g. y* v' _
p=0.0;
3 N1 t* ~ X/ w* y* S3 Z
for (i=0;i<=n-1;i++)
6 Y' s/ v; F5 ]/ G
{
8 S) P$ J% n) q% G8 w2 O
x=a+(i+0.5)*h;
, n ]- Y# r6 j- S0 i$ A# c
p=p+Func(x);
7 n$ u5 U0 E D1 v3 u
}
+ h. h8 S; G$ q; d) j" H% a/ ^( K; z
8 U5 @" x; e8 @1 P% v' q
p=(y[0]+h*p)/2.0;
- F: L4 G* W: k* n
s=1.0;
6 g, ^( v1 M- {1 C$ g, ]. C
for (k=1;k<=m;k++)
. n) E: B5 W" B1 t
{
, [8 c, {+ p! b; ?; A2 ?
s=4.0*s;
! A: d7 k8 y/ K5 k- ~
q=(s*p-y[k-1])/(s-1.0);
! }4 A' r" [- g% G; \' J
y[k-1]=p; p=q;
+ D2 z' F1 r8 X+ u& u- ?) O% R
}</P>
: L& i5 d2 J7 ]3 R; H, W
<
> ep=fabs(q-y[m-1]);
5 x3 ?/ d. A( o/ y8 A; ]8 W" c: {
m=m+1;
# g, @! b/ ~+ d* m, U
y[m-1]=q;
/ _" {8 z! |2 U
n=n+n;
) S- b* T5 N: i7 ~/ N9 U
h=h/2.0;
/ \- [2 S$ Y8 a% w1 T7 F
}
, `% \0 ~8 ^& P& B% j( W
L2 m8 r! c9 p* Z1 S3 F
return(q);
* D, m( R- Q/ S+ p8 R9 `
}</P>
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
>积分的源代码,效率很好的> // 迭代初值> ep=fabs(q-y[m-1]);5 x3 ?/ d. A( o/ y8 A; ]8 W" c: {