实验二、定积分的近似计算

将自变量看成是向量
20
quad 举例
例：用 quad 计算定积分： I
1 dx 0 1 x2
解： >> quad('1./(1+x.^2)',0,1)
>> quad('1./(1+x.^2)',0,1,1e-10)
>> quad('1./(1+x.^2)',0,1,1e-16)
函数表达式一定要用 单引号 括起来！ 涉及的运算一定要用 数组运算！
计算函数 f 关于自变量 v 的不定积分
findsym(f,1
int(f,a,b)
)
计算 f 关于默认自变量 的定积分，积分区间为[a,b]。
int(f) 计算 f 关于默认自变量 的不定积分。
25
int 举例
例：用 int 函数计算定积分：
I
1 dx 0 1 x2
解： >> syms x; >> f=1/(1+x^2); >> int(f,x,0,1)
中点法：
1 dx 0 1 x2
n
h
i 1
f ( xi1 2
xi )
0.78540024673078
6
矩形法举例
相对误差分析
理论值：
1 dx 0 1 x2

arctan x
1 0

π 4
左点法相对误差：0.78789399673078 / 4 0.003178 /4
例：计算二重积分 I 2 1 (4xy 3x2 )dxdy 0 1 >> dblquad(inline('4*x*y+3*x.^2') ,-1, 1, 0, 2) >> dblquad(@(x,y)4*x*y+3*x.^2 , -1, 1, 0, 2)
指定 x、y 分别是第一和第二积分变量

y2n
4( y1 2( y2

y3 y4

0.78539816339745
y2n1 ) y2n2 )]
2.827 10 相对误差：0.78539816339745 / 4
-16
/4
16
Matlab 函数
数值积分的常见算法
矩形法 梯形法 抛物线法
(y0 4y1 y2 )
ba 6n
(y0 4 y1 y2 )
13
抛物线法
同理可得：
x4 x2
f ( x)dx

ba 6n
(
y2
4 y3

y4 )
x2n x2 n 2
f ( x)dx

ba 6n ( y2n2
4 y2n1

y2n )
相加即得：
抛物线法公式 或 辛卜生 (Simpson) 公式
15
抛物线法
例：用抛物线法计算下面定积分 ( 取 n=100 )，并计算相对误差
I 1 dx 0 1 x2
解： a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 )
==>
1 dx 0 1 x2

ba 6n [ y0
h ba n
点 i [ xi1, xi ] 可以任意选取，常见的取法有： 左端点 xi1 , 右端点 xi 和中点( xi1 xi ) / 2 。
左点法 右点法
中点法
4
矩形法
b
n
f ( x)dx
a
f (i )xi
i 1
x1 x2
xi
xn
x0 步长
n
Si
i 1

n i 1
yi1 2
yi xi
Si
9
梯形法
如果我们 n 等分区间 [a,b]，即令：
x1 x2 xn
h ba n
则 S
b a
f ( x)dx
n
Si
i 1

n i 1
yi1 2
yi xi
n
h
i 1
yi1 2
例：计算二重积分 I 2 1 (4xy 3x2 )dxdy 0 1 X >> dblquad(inline('4*x*y+3*x^2'),-1,1,0,2) >> dblquad(inline('4*x*y+3*x.^2'),-1,1,0,2) 23
dblquad 举例
被积函数 f (x,y) 的另一种定义方法：匿名函数
11
抛物线法
2n 等分区间 [a,b] ，得
h1

ba, 2n
xi

ih1 ,
i

0, 1,
, 2n
计算每个节点上的函数值：
yi f ( xi ), i 0,1, , 2n
在区间 [x0, x2] 上，用过以下三点 P0 ( x0 , y0 ), P1( x1, y1 ), P2( x2, y2 )
定积分的近似计算
问题背景和实验目的
定积分计算的基本公式是牛顿－莱布尼兹公式。但当 被积函数的原函数不知道时，如何计算？这时就需要利 用近似计算。特别是在许多实际应用中，被积函数甚至 没有解析表达式，而是一条实验记录曲线，或一组离散 的采样值，此时只能用近似方法计算定积分。
本实验主要研究定积分的三种近似计算算法：矩形法、 梯形法和抛物线法。同时介绍 Matlab 计算定积分的相 关函数。
i 1
i 1
右点法
b a
f ( x)dx
n i 1
f
(
xi1 2
xi
)xi
n
h
i 1
f ( xi1 xi ) 2
中点法
5
矩形法举例
例：用不同的矩形法计算下面的定积分 ( 取 n=100 )，
并比较这三种方法的相对误差。
I
1 dx 0 1 x2
解：a=0, b=1, n=100
的抛物线来近似原函数 f (x) 。
用抛物线代替该直线， 计算精度是否会更好？
12
抛物线法
设过以上三点的抛物线方程为：
y = x2 + x + = p1(x)
则在区间 [x0, x2] 上，有
x2 f (x)dx x0
x2 x0
p1( x)dx
x2 ( x2 x )dx
1
主要内容
数值积分的常见算法
矩形法 梯形法 抛物线法
Matlab 求积分函数
数值积分函数：trapz、quad、dblquad 符号积分函数：int
2
定积分的近似
定积分的定义
b
f ( x)dx
a
n
lim
n
f (i )xi , i [xi1, xi ]
按字母顺序，大写字母排在小写字母的前面
22
dblquad 举例
f 中关于第一自变量的运算是数组运算，即把 x 看成是向 量，y 看成是标量。也可以全部采用数组运算
例：计算二重积分 I 2 1 (4xy 3 y2 )dxdy 0 1 >> f=inline('4*x*y+3*y^2'); >> I=dblquad(f,-1,1,0,2)
b
n
f ( x)dx
x2i f ( x)dx
a
i1 x2i2

n i 1
ba 6n
( y2i2
4 y2i1

y2i )
14
抛物线法
整理后可得：
b
f ( x)dx
a

ba 6n
[
y0

y2n

4(
y1

y3

2( y2 y4
y2n1 ) y2n2 )]
x1
x2
xi

h

b
n
a
xi1 xi
xi a ih, i 1,2,
xn1
xn
n
节点
b
n
n
f ( x)dx
a
f ( xi-1 )xi h f ( xi1 )
i 1
i1
左点法
b
n
n
f ( x)dx
a
f ( xi )xi h f ( xi )
h =1/n=0.01, xi = i*h,
(i = 0, 1, 2, ..., 100)
左点法：
1 dx 0 1 x2
n
h
i 1
f ( xi1 )
0.78789399673078

右点法：
1 0
1
dx x2
n
h
i 1
f ( xi )
0.78289399673078
==> h =1/100=0.01, xi = i*h, yi = f (xi)
==>
1 dx 0 1 x2

h

y0 2

y1


yn1

yn 2

0.78539399673078
相对误差： 0.78539399673078 / 4 5.305 10-6 /4
右点法相对误差：0.78289399673078 / 4 0.003188 /4
中点法相对误差：0.78540024673078 / 4 2.653 10-6 /4
不同的算法有不同的计算精度
有没有更好的近似计算定积分的方法 ?
7
定积分几何意义
b
S a f ( x)dx
下面的命令运行结果和上面的一样吗？
>> dblquad(@(y,x)4*x*y+3*x.^2 , -1, 1 , 0, 2 )
24
int
b
a f (v)dv
f (v)dv
符号积分： int
int(f,v,a,b) 计算函数 f 关于自变量 v 的定积分，积分区间为 [a, b]
int(f,v)
x0 i1
x1 x2
x0
x1
x2
xi xi1 xi
xn xn1
xn
xi xi xi1
x

max i
xi
3
矩形法
定积分的近似：
b
n
f ( x)dx
a
f (i )xi
i 1
n 充分大，x 充分小
通常我们取 x1 x2 xn
a f ( x)dx
quad(f,a,b,tol)
f = f(x) 为被积函数，[a,b] 为积分区间，tol 为计算精度
不用自己分割积分区间
可以指定计算精度，若不指定，缺省精度是 10-6
精度越高，函数运行的时间越长
此处的函数 f 是数值形式，应该使用数组运算，即： .* ./ .\ .^
前面的 做法
1 dx 0 1 x2

ba n
y0 2

y1

y2


yn1
yn 2

trapz 函数
>> x=0:1/100:1; >> y=1./(1+x.^2); >> trapz(x, y)
trapz(x,1./(1+x.^2))
19
quad
ቤተ መጻሕፍቲ ባይዱuad
抛物线法
b
yi
==>
b a
f
( x)dx

h
y0 2

y1


yn1

yn 2

梯形公式
梯形公式与中点公式有什么区别 ?
10
梯形法举例
例：用梯形法计算下面定积分 ( 取 n=100 )，并计算相对误差
I
1 dx 0 1 x2
解： a=0, b=1, n=100, f (x) = 1/( 1+x2 )
Matlab 求积分函数
数值积分函数：trapz、quad、dblquad 符号积分函数：int
17
trapz
trapz
梯形法
trapz(x,y)
x 为分割点（节点）组成的向量， y 为被积函数在节点上的函数值组成的向量。
b a
f
( x)dx

b
n
a

y0 2

y1


y
f (x)

S1 S2
Si
Sn
oa
xi1 xi
b
b
n
S f ( x)dx a
Si
i 1
x
8
梯形法
曲边小梯形的面积可以由直边小梯形的面积来近似
整个曲边S梯i 形y的i1面2积yi xi yi f ( xi ), i 1, 2, , n
：
b
S a f ( x)dx
x0

x3
3 x2
6


x2


x
x2

2
x0
(x02
x0
x0
3 )
(
x23

x03
)

2
(
(x22 x2
x22 )

x02
)

( x2 x0 )2 2( x2 x0
( x2 ) 4
x0
)
x2 x0 6
21
dblquad
db
c a f ( x, y)dxdy
抛物线法计算二重积分： dblquad
dblquad(f,a,b,c,d,tol)
tol 为计算精度，若不指定，则缺省精度为 10-6
f 可以是： 字符串；inline 定义的内联函数；函数句柄
[a,b] 是 第一积分变量 的积分区间， [c,d] 是 第二积分变量 的积分区间
yn1

yn 2

x [ x0 , x1, , xn] y [ f ( x0 ), f ( x1), , f ( xn )]
18
trapz 举例
例：用梯形法计算下面定积分 ( 取 n=100 )
I 1 dx 0 1 x2
解： a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 )