数学实验题目2 Romberg积分法
巡山小妖精
779次浏览
2021年01月30日 01:13
最佳经验
本文由作者推荐
挫折的作文-圣诞节的祝福语
数学实验题目
2 Romberg
积分法
摘要
考虑积分
I
(
f)
欲求其近似值，可以采用如下公式：
n
1
b
a
f
(
x
)
dx
h< br>[
f
(
x
i
)
f
(
x< br>i
1
)]
i
0
2
b
a
2
h
f
(
)
[
a
,
b
]
E
12
n
1
h
（复化）辛卜生公式
S
[
f
(
x
i
)
4f
(
x
1
)
f
(
x
i
1
)]
i
i
0
6
2
（复化）梯形公式
T
b
a
h
(4)
E
f
(
)
[
a
,
b
]
180
2
n
1
h
（复化）柯特斯公式
C
[7
f
(
x
i
)
32
f
(
x
1
)
12
f
(
x
1
)
i
i
i
0
90
4
2
32
f
(
x
i
3
4
4
)
7
f
(
x
i
1
)]
6
2(
b
a
)
h
(6)
E
f
(
)
[
a
,
b
]
945
4
这里，梯形公式显得算法简单，具有如下递推关系
1h
n
1
T
2
n
T
n
f
(
x
1
)
i
2
2
i
0
2
因此，
很容易实现从低阶的计算 结果推算出高阶的近似值，
而只需要花费较少的附加函数计
算。但是，由于梯形公式收敛阶较低 ，收敛速度缓慢。所以，如何提高收敛速度，自然是人
们极为关心的课题。为此，记
T
0,
k
为将区间
[
a
,
b
]
进行
2
等份的复化梯形积分结果，
T
1,
k
为
将区间
[
a
,
b
]
进行
2
等份的复化辛卜生积分结果，T
2,
k
为将区间
[
a
,
b
]
进行
2
等份的复化柯
特斯积分结果。根据李查逊（
Richardson< br>）外推加速方法，可得到
k
k
k
T
m
,< br>k
4
m
T
m
1,
k
1
T
m
1,
k
4
m
1
k
0,1,2,
m
0,1,2,
可以证明，如果
f
(
x
)
充分光滑，则有
lim
T
m
,
k
I
(
f
)
（
m
固定）
k
lim
T
m
,0
I
(f
)
m
这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
该方法的计算可按下表进行
T
0,0
T
0
,
1
T
0
,
2
…
T
0,
m
T
1,0
T
1
,
1,
m
1
1
…
T
T
2,0
…
T
2,
m
2
…
…
T
m
,0
很明显，龙贝格计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果
T
m< br>,
k
，可存放在
T
0,
k
位置上，其最终结果
T
m
,0
是存放在
T
0,0
位置上。
前言
利用龙贝格积分法计算积分
I
(
f
)
b
a
f
(
x
)
dx
程序设计流程：
1
．准备初值，计算
T
0,0
且
k< br>
0
（
k
为等份次数）
2
．按梯形公式的递推关系，计算
a
b[
f
(
a
)
f
(
b
)]< br>
2
1
b
a
2
1
b< br>
a
1
T
0,
k
1
T
0,
k
k
1
f
(
a
k
(
i
))
2
2
2
2
i
0
3
．按龙贝格公式计算加速值
k
T
0,
k
m
T
m
,
k
m
4
m
T
m
1,
k
1
m
T
m
1,
k
m< br>4
1
m
m
0,1,2,
,< br>k
4
．精度控制。对给定的精度
，若
|
T
m
,0
T
m
1,0
|
则终止计算，并取
T
0,
s
T
m,
s
作为所求结果；否则
k
k
1
，重复
2~4
步，直到满足精度
为止。
问题
1
（
1
）
程序运行如下：
I = Romberginterg(inline('x.^2.*exp(x)'),0,1,25,1e-6)
I = 0.7183
（
2
）
程序运行如下：
I = Romberginterg(inline('exp(x).*sin(x)'),1,3 ,25,1e-6)
I = 10.9502
（
3
）
程序运行如下：
I = Romberginterg(inline('4./(1+x.^2)'),0,1,25,1e-6)
I = 3.1416
（
4
）
程序运行如下：
I = Romberginterg(inline('1./(1+x)'),0,1,25,1e-6)
I = 0.6931
问题
2
（
1
）
程序运行如下：
I = Romberginterg(inline('cos( x)./sqrt(1-x)'),0,1,25,1e-6)
I = NaN
因为函数在
x=0
处出现了
0/
0
的情况，极限为
1
，所以< br>Matlab
的结果显示为
NaN
非数，
解决方法是把下限
0
改为一个小数
eps
。
I = Romberginterg(inline('sin(x)./x'),eps,1,25,1e-6)
I = 0.9461
（
2
）
程序运行如下：
I = Romberginterg(inline('cos(x)./sqrt(1-x)'), 0,1,25,1e-6)
I = NaN
与
（
1
）
的 原因相同，
函数在
x=1
处出现了
0/
0
的情况，
结果为无穷，
此时应选择变换
，
将积分变为
，再进行变换
，将积分变 为
，变换后输入命令：
I = RombergInterg(inli ne('2*sin(x).*cos((sin(x)).^2)'),0,pi/
2,25,1e- 6)
I = 1.4996
（
3
）
程序运行如下：
I = Romberginterg(inline('cos( x)./sqrt(x)'),0,1,25,1e-6)
I = NaN
函数在
x=0
处出现了
1
/
0
的情况，结果为无穷。先将积分变为
，再做
变换
，
I = 2*Romberginterg(inline('cos(x.^2)'),0,1,25,1e-6)
I = 1.8090
（
4
）
程序运行如下：
I = Romberginterg(inline('x.*sin(x)./(1-x.^2)' ),-1,1,25,1e-6)
I = NaN
函数在
x=1,-1
处 出现了
sin1
/
0
的情况，
结果为无穷。
被积函数为偶函 数，
做变换
，
积分变为
I = 2*Romberginterg (inline('sin(x).*sin(sin(x))'),0,pi/
2,25,1e-6)
I = 1.3825
或者使用
Gauss- Chebyshev
求积公式：
I = GaussInterg(inline( 'x.*sin(x)'),'Chebyshev',-1,1,1e-6)
I = 1.3825
由于
Gauss-Chebyshev
求积公式求的是
。
在
的积分因此此处的
为
问题
3
（
1
）
程序运行如下：
I = Romberginterg(inline('exp(-x.^2)'),-10,10,25,1e-6)
I = 1.7725
由于积分区间无限，而函数在
[-10,10]
外很 小，可以用函数在
[-10,10]
上的积分近似这个
无穷积分。或者使用
G auss-Hermite
求积公式：
I = GaussInterg(inline('x+1-x'),'Hermite',0,0,1e-6)
I = 1.7725
由于
Gauss- Hermite
求积公式求的是
（
2
）
程序运行如下：
先做变换
，则原积分变为
在
的积分因此此处的
为
1
。
I = 2*Romberginterg(inline('1./(1+x.^2)'),0,1,25,1e-6)
I = 1.5708
（
3
）
程序运行如下：
使用
Gauss-Hermite
求积公式：
I = Gauss Interg(inline('cos(x).^3'),'Hermite',0,0,1e-6)
I = 1.0820
由于
Gauss- Hermite
求积公式求的是
在
的积分因此此处的
为
。
（
4
）
程序运行如下：
使用
Gauss-Laguerre
求积公式：
I = Gaus sInterg(inline('sin(x).^2'),'Laguerre',0,0,1e-6)
I = 0.4000
由于
Gauss- Laguerre
求积公式求的是
在
的积分因此此处的
为
。
使用的函数
function I = RombergInterg(fun, a, b, npanel, tol, flag)
% RombergInterg
用
Romberg
方法求积分
%
% Synopsis:
I = RombergInteg(fun,a,b,n,tol)
%
% Input:
fun
= (string)
被积函数的函数名
%
a, b
=
积分下限和积分上限
%
npanel = (optional)
将积分区间平分的段数，默认为
25
%
tol
= (optional)
计算误差上限，默认为
5e-9
%
flag
= (optional)
是否显示
Romberg-T
数表，默认为
0
为不显示
%
% Output:
I =
通过
Romberg
方法求积分的近似值
if nargin < 4
npanel = 25;
end
if nargin < 5
tol = 5e-9;
end
if nargin < 6
flag = 0;
end
T(1,1) = TrapezoidInteg(fun, a, b, npanel);
%T0(h) = T(h)
err = 1;
%
初始化误差值
m = 2;
%
初始化行计算的行序号
while err >= tol
T(m,1) = TrapezoidInteg(fun, a, b, 2^m*npanel);
%
计算第
m
行
T0(h/
2^m)
T(m,m) = 0;
%
将矩阵
T
变成
m*m
for n = 2:m
T(m,n)
=
(
4^(n-1)*T(m,n-1)
-
T(m-1,n-1))
/
(4^(n-1)
-
1);
%Tm(h/
2^k)
与
Tm-1(h/
2^(k+1) )
和
Tm-1(h/
2^k)
的递推关系
end
err = abs( T(m,m) - T(m-1,m-1) );
%
计算误差值
m = m + 1;
%
计算下一行
end
I = T(m-1,m-1);
if flag ~= 0
disp(T);
end
function
I = TrapezoidInteg(fun, a, b, npanel)
% TrapezoidInteg
用复化梯形公式求积分
%
% Synopsis:
I = TrapezoidInteg(fun,a,b,n)
%
% Input:
fun
= (string)
被积函数的函数名
%
a, b
=
积分下限和积分上限
%
npanel = (optional)
将积分区间平分的段数，默认为
25
%
% Output:
I =
通过复化梯形公式求积分的近似值
if nargin < 4
npanel = 25;
end
nnode = npanel + 1;
%
节点数
=
段数
+ 1
h = (b-a)/(nnode-1);
%
步长
x = a:h:b;
%
将积分区间分段
f = feval(fun,x);
%
求节点处被积函数的值
I = h * ( 0.5*f(1) + sum(f(2:nnode-1)) + 0.5*f(nnode) );
function I = GaussInterg(fun, type, a, b, tol)
% GaussInterg
用
Gauss
型求积公式求积分
,
具体形式由使用者选取
%
% Synopsis:
I = GaussInterg(fun, type, a, b)
%
I = GaussInterg(fun, type, a, b, tol)
%
% Input:
fun
= (string)
被积函数的函数名
%
type = (string)
具体
Gauss
求积公式的形式
%
a,b
=
积分上下限，
Laguerre
只计算0
到
inf
，
Hermite
只计算
-inf
到
inf
，所
以
a,b
对这两种形式无效
%
tol
= (optional)
误差容忍限度，默认为
5e-5
%
% Output:
I =
通过
Gauss
型求积公式求积分的近似值
if nargin < 5
tol = 5e-5;
end
n = 7;
%
默认从
7
节点多项式开始计算
IOld = 1;
%
初始化
IOld
为
1
err = 1;
%
初始化误差为
1
while err >= tol
switch type
case 'Legendre'
%
计算无奇点被积函数在
-1
到
1
的积分
I = GaussLegendreInterg(fun, a, b, n);
case 'Chebyshev'
%
计算被积函数形如
f(x)/sqrt(1-x^2)
在
-1
到
1
的积分
I = GaussChebyshevInterg(fun, a, b, n);
case 'Laguerre'
%
计算被积函数形如
exp(-x)*f(x)
的在
0
到
inf
的积分
I = GaussLaguerreInterg(fun, n);
case 'Hermite'
%
计算被积函数形如
exp(-x^2)*f(x)
的在
-inf
到
inf
的积分
I = GaussHermiteInterg(fun, n);
otherwise
error('No such type!');
end
err = abs(I - IOld);
%
计算误差
IOld = I;
%
把
IOld
赋值为
I
进行下次迭代
n = n+1;
%
多项式节点递增
end