Posted on 2006-05-29 09:25
Tauruser 阅读(1768)
评论(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下编译运行通过。