验证数值积分求积公式及复合梯形公式程序设计

《复合梯形公式》实验报告
实验名称：验证数值积分求积公式及复合梯形公式程序设计成绩：___________
专业班级：数学与应用数学1202班姓名：王晓阳学号：2012254010228 实验日期：2014 年10月20日
实验报告日期：2014年11月3日
一．实验目的
1掌握定积分的数值求解方法，验证数值积分求积公式.
2.掌握数值积分的基本思想.
3.掌握matlab实现数值积分函数的调用格式.
4.编写复合梯形公式matlab程序及学会调用.
5.学会用复合梯形公式求函数近似解.
二、实验内容
1.数值积分的实现
(1)被积函数是一个解析式
Matlab提供了quad函数和quadl函数来求定积分.它们的调用格式为：
Quad(filename,a,b,tol,trace)
Quadl(filename,a,b,tol,trace)
其中filename 是被积函数名。

a 和b 分别是定积分的下限和上限。

Tol 用来控制积分精度，默认时取610tol -=。

Trace 控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0。

例6.20 用两种不同的方法求2
10x I e dx -=⎰. (2)被积函数由一个表格定义
在matlab 中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。

其中向量X 、Y 定义函数关系Y=f(X).X 、Y 是两个等长的向量；
()12,,n X x x x =，()12,,n Y y y y =，并且12n x x x <<<，积分区间是[]1,n x x 。

例6.21用trapz 函数计算210
x I e dx -=⎰. (3)二重积分数值求解
Matlab 提供的dblquad 函数可以直接求出二重定积分的数值解。

该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)
该函数求f(x,y)在[a,b]*[c,d]区域上的二重定积分。

例6.22 计算二重定积分()2
122212sin x e x y dxdy ---+⎰
⎰
2.复合梯形公式 由于牛顿-柯特斯公式在8n ≥时不具有稳定性，故不可能通过提高阶的方法来提高求积精度.为了提高精度通常可把积分区间分成若干子区间(通常是等分)，再在每个子区间上用低阶求积公式，这种方法称为复合求积法.
将区间[],a b 划分为n 等份，分点0,0,1,
k x x kh k n =+=，0,,n b a a x b x h n
-===
，在每个子区间[]1,k k x x +上采用梯形公式
()()()2b
a b a f x dx f a f b -≈
+⎡⎤⎣⎦⎰，则得()()()()()1
111002k k n n b x k k n a x k k h I f x dx f x dx f x f x R f +--+=====++⎡⎤⎣
⎦∑∑⎰⎰. 记()()()()()11101222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎡⎤⎣⎦⎢⎥⎣⎦
∑∑，称为复合梯形公式，其余项可由[]
()()()3'',,12b a R f f a b ηη-=-∈得[]()()3
1''10,,12n n n k k k k k h R f I T f x x ηη-+=⎡⎤=-=-∈⎢⎥⎣⎦
∑.由于
()}{()()}{
1''''''01min max n k k k k f f f n ηηη-=≤≤∑，所以(),a b η∃∈，使()()1''
''01n k k f f n ηη-==∑.于是复合梯形公式的余项为[]()2''12
n b a R f h f η-=-.(),a b η∈. ()lim 0n n R f →∞=即lim n n T I →∞
=，所以复合梯形公式收敛. 以以原函数()()sin x f x x
=为例，进行matlab 程序设计. 三、实验环境
MATLAB R2014a
四、实验步骤
(一)验证数值积分公式
(1)被积函数是一个解析式
例6.20先建立一个函数M文件ex.m:
function ex=ex(x)
ex=exp(-x.^2);
end
然后在matlab命令窗口，输入命令：
format long
I=quad('ex',0,1)
(2)被积函数由一个表格定义
例6.21
在matlab命令窗口，输入命令：
X=0:0.01:1;
Y=exp(-X.^2);
trapz(X,Y)
(3)二重积分数值求解
例6.22
(1)建立一个函数文件fxy.m
function f=fxy(x,y)
global ki;
ki=ki+1; %ki用于统计被积函数的调用次数f=exp(-x.^2/2).*sin(x.^2+y);
end
(2)调用dblquad函数求解
global ki;ki=0;
I=dblquad('fxy',-2,2,-1,1)
(二)复合梯形公式算法流程图
(三)程序设计
以原函数()
() sin x
f x
x
=为例，a,b分别为积分区间的上限和下限，n为将积分区间划分的份数.
1．编写原函数M文件f.m
function y=f(x)
y=sin(x)/x;
2.编写求复合梯形公式M文件FT.m
function Tn=FT(a,b,n)
a=input('请输入a的值：');
b=input('请输入b的值:');
n=input('请输入n的值:');
h=(b-a)/n;
for k=0:n
x(k+1)=a+k*h;
if x(k+1)==0
x(k+1)=10^(-10);
end
end
T1=h/2*(f(x(1))+f(x(n+1)));
for i=2:n
F(i)=h*f(x(i));
end
T2=sum(F);
Tn=T1+T2;
3.运行程序
4.用quad求解与复合梯形公式求解结果做对比I=quad('sin(x)./x',0,1)
五、实验结果
(一) 验证数值积分公式
(1) 被积函数是一个解析式
I =0.746824180726425
(2) 被积函数由一个表格定义
ans =0.7468
(3) 二重积分数值求解
I =1.5745
(二)复合梯形公式程序运行结果：
FT
请输入a 的值：0
请输入b 的值：1
请输入n 的值：8
ans =0.9457
FT
请输入a 的值：0
请输入b 的值：1
请输入n 的值：3
ans =0.9433
调用quad 命令运行结果
I =0.9461
六、实验讨论、结论
1.对于函数()()sin x f x x
=，将积分区间[0,1]划分为8等份，应用复合梯形公式matlab 程序求得积分()10sin x I dx x
=⎰的近似值80.9457T =.调用quad 命令求得I=0.9461，余项[]n n R f I T =-为0.0004.当将积分区间划分为3等份时，30.9433T =，余项为0.0028，由此可知，将区间划分的份数越多，结果越准确.
2.通过实验进一步学习了积分的数值计算方法，设法构造某个简单函数近似
表示，然后得到近似值。

公式不但具有最高代数精度，而且收敛性和稳定性都有保证，因此是高精度的求积公式。

七、参考资料
[1]李庆杨，王能超，易大义. 数值分析. 清华大学出版社，2008，P106
[2]刘卫国. MATLAB程序设计与应用. 高等教育出版社，2006，P157。