Posted on 2006-05-29 09:25
Tauruser 閱讀(1807)
評論(1) 編輯 收藏 引用 所屬分類:
數值計算
先介紹一下Romberg求積。
6.3 外推原理與Romberg求積
6.3.1 復合梯形公式遞推化與節點加密
在計算機上用等距節點求積公式時，若精度不夠可以逐步加密節點.設將區間
分為n等分，節點
，在區間
上梯形公式為
若節點加密一倍，區間長為
，記中點為
在同一區間上的復合梯形公式惟
于是
(6.3.1)
它表明
是在
的基礎上再加新節點
的函數值之和乘新區間長
，而不必用(6.2.6)重新計算，這時有誤差估計式
若
，則得
(6.3.2)
它表明用
，其誤差近似
.這也是在計算機上估計梯形公式誤差的近似表達式.若
(給定精度)，則
.
若在區間［a,b］中做2n等分時，在上用Simpson公式計算，則由(6.2.8)可知
它恰好是(6.3.2)中I(f)的近似值，即
它表明用(6.3.2)計算I(f)，其精度已由
提高到
如果再將區間分半，使分為4個小區間，長度為
，則可由(6.3.1)計算出
及
，利用復合公式余項(6.2.9)得
如果
，則有
(6.3.3)
從而有復合Simpson公式的誤差估計
如果用(6.3.3)近似
，即
(6.3.4)
則精度可達到
.類似做法還可繼續下去.這樣對區間逐次分半，利用公式(6.3.1)逐次遞推.再由(6.3.2)，(6.3.3)逐次構造出精度愈來愈高的計算積分I(f)的公式，這就是Romberg求積的基本思想.
以下為我自己寫的求積程序。
//
?RombergIntegral.cpp?:?定義控制臺應用程序的入口點。
//
#include?
<
cmath
>
#include?
<
iostream
>
#include?
<
vector
>
using
?
namespace
?std;
const
?
double
?PRECISION(.
000001
);
//
精度控制
const
?unsigned?
int
?MAXK(
20
);
//
求解步驟控制
double
?RombergIntegral(
double
?(
*
f)(
double
?x),
double
?a,?
double
?b);
vector
<
vector
<
double
>>
?T;
//
用于存儲T表
double
?f(
double
?x)
//
要求的積分函數
{
????
return
?x
*
sin(x);
}
int
?_tmain(
int
?argc,?_TCHAR
*
?argv[])
{
????cout
<<
"
本程序用于求解函數f(x)=x*sin(x)在0到6.28的積分
"
<<
endl;
????cout
<<
"
積分結果為：
"
<<
RombergIntegral(f,
0
,
6.28
)
<<
endl;
????cout
<<
"
精度為
"
<<
PRECISION
<<
endl;
????
return
?
0
;
}
double
?RombergIntegral(
double
?(
*
f)(
double
?x),
double
?a,?
double
?b)
{
????
int
?k(
0
);
????
double
?h
=
b
-
a;
????vector
<
double
>
?temp;
????T.push_back(temp);
????T[
0
].push_back(h
*
((
*
f)(a)
+
(
*
f)(b))
/
2
);
????
for
(k
=
1
;
1
;
++
k)
????
{
????????T.push_back(temp);
????????T[
0
].push_back(
0.5
*
T[
0
][k
-
1
]);
????????
for
(
int
?i
=
0
;i
<
pow(
2
.,k
-
1
);
++
i)
????????
{
????????????T[
0
][k]
+=
0.5
*
h
*
((
*
f)(a
+
h
/
2
+
i
*
h));
????????}
????????
for
(
int
?i
=
1
;i
<=
k;
++
i)
????????????T[i].push_back((pow(
4
.,i)
*
T[i
-
1
].back()
-
T[i
-
1
][T[i
-
1
].size()
-
2
])
/
(pow(
4
.,i)
-
1
));
????????h
/=
2
;
????????
double
?temp
=
T[k].back();
????????
if
(fabs(T[k].front()
-
T[k
-
1
].front())
<
PRECISION?
||
??k
==
MAXK)?
break
;
//
到
????}
????
????
return
?T[k].back();
}
//
以上程序在vs2005+win2003下編譯運行通過。