Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
积分
Presentation
?首页：题目，作者，时间，单位，资助
?介绍：以前工作的介绍，背景，经验和教训，参考文献
?介绍二：你的工作和主要结果，和比较
?主要内容：理论、计算、试验，困难，有趣的事情
?总结：主要的结论，意义，不足
?展望：可能的发展，前景，应用，可以深入的工作
用图形、动画、多媒体，美工
详尽的资料，观众的预期
0 1 2 3 4 5 6
0
2
4
6
8
10
12
x
x (2+sin(2 x))
Quadrature
?ba dxxf )(
Quadrature
??? ?? bccaba dxxfdxxfdxxf )()()(
基本求积公式
?????? ???? 2 bahfMfd xba
中点公式：
梯形公式 (trapezoid rule)
2
)()( bfafhTf d xb
a
????
The order of a quadrature rule is the degree of the lowest
degree polynomial that the rule does not integrate exactly.
两阶求
积公式
0 1 2 3 4 5 6
0
2
4
6
8
10
12
x
x (2+sin(2 x))
Adaptive Quadrature
? ?10 2 31dxx
4
1
2
11 2 ??
?
??
?
??M
2
1
2
101 2 ???T
12
1
4
1
3
1 ??
6
1
2
1
3
1 ???
)(2)( MSTS ????
TMS 3132 ??
Simpson’s rule:
2) ),()(4)((6
bacbfcfafhS ?????
TMS 3132 ??
四阶求积公式
Composite quadrature rule,对 [a,c],[c,d],d=(a-c)/2,e=(c-b)/2
)),()(4)(2)(4)((122 bfefcfdfafhS ?????
Adaptive Quadrature
S和 S2都逼近同样的积分，他们的差可以用来估计误差
)( 2 SSE ??
二者也可以组合起来得到更精确的逼近。这两个求积
公式都是四阶的，但 S2的步长是 S的一半
)(16)( 2SQSQ ???
15/)( 22 SSSQ ???
六阶求积公式
Weddle’s rule,the sixth order Newton-Cotes rule or the
first step of Romberg integration.
Extrapolated Simpson’s rule
quadtx
Matlab中 quad用递归自适应的外插 Simpon公式
quadtx是一个简化版本。先求出三点的函数值，然后用递
归过程 quadtxstep完成积分。
quadtx.m
[Qa,ka] = quadtxstep(F,a,c,tol,fa,fd,fc,varargin{:});
[Qa,ka] = quadtxstep(F,c,b,tol,fc,fe,fb,varargin{:});
tol的故事：
http:// www.inf.ethz.ch/personal/gander,
Specifying Integrands
? ?10 41 1 dxx
f = inline(‘1/sqrt(1+x^4)’)
Q = quadtx(f,0,1)
?10 )s i n ( dxx x
f = inline(‘sin(x)/x’)
Q = quadtx(f,0,pi)
Q = quadtx(f,realmin,pi)
Q = quadtx(@sinc,0,pi)
? ?? ?10 11 )1( dxtt z ? f = inline(‘t^(z-1)*(1-t)^(w-1)’,’t’,’z’,’w’)
z = 8/3; w = 10/3; tol = 1.e-6;
Q = quadtx(f,0,1,tol,z,w)
vectorize
10
20
30
40
50
60
70
80
90
100
93
Performance
04.0)9.0(
1
01.0)3.0(
1)(
22 ?????? xxxh
humps 函数
quadgui(@humps,0,1,1.e-4)
Performance
syms x
h = 1/((x-.3)^2+.01) + 1/((x-.9)^2+.04) - 6
I = int(h)
D = simple(int(h,0,1))
Qexact = double(D)
for k=1:12
tol=10^(-k); [Q,fcount]=quadtx(@humps,0,1,tol);
err=Q-Qexact; ratio=err/tol;
fprintf('%8.0e %21.14f %7d %13.3e %9.3f\n',...
tol,Q,fcount,err,ratio),end
tol Q count err ratio
1e-001 29.83328444174864 25 -2.504e-002 -0.250
1e-002 29.85791444629948 41 -4.109e-004 -0.041
1e-003 29.85834299237637 69 1.760e-005 0.018
1e-004 29.85832444437543 93 -9.511e-007 -0.010
1e-005 29.85832551548643 149 1.200e-007 0.012
1e-006 29.85832540194041 265 6.442e-009 0.006
1e-007 29.85832539499819 369 -5.005e-010 -0.005
1e-008 29.85832539552631 605 2.763e-011 0.003
1e-009 29.85832539549604 1061 -2.640e-012 -0.003
1e-010 29.85832539549890 1469 2.274e-013 0.002
1e-011 29.85832539549867 2429 -7.105e-015 -0.001
1e-012 29.85832539549867 4245 0.000e+000 0.000
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
积分
Presentation
?首页：题目，作者，时间，单位，资助
?介绍：以前工作的介绍，背景，经验和教训，参考文献
?介绍二：你的工作和主要结果，和比较
?主要内容：理论、计算、试验，困难，有趣的事情
?总结：主要的结论，意义，不足
?展望：可能的发展，前景，应用，可以深入的工作
用图形、动画、多媒体，美工
详尽的资料，观众的预期
0 1 2 3 4 5 6
0
2
4
6
8
10
12
x
x (2+sin(2 x))
Quadrature
?ba dxxf )(
Quadrature
??? ?? bccaba dxxfdxxfdxxf )()()(
基本求积公式
?????? ???? 2 bahfMfd xba
中点公式：
梯形公式 (trapezoid rule)
2
)()( bfafhTf d xb
a
????
The order of a quadrature rule is the degree of the lowest
degree polynomial that the rule does not integrate exactly.
两阶求
积公式
0 1 2 3 4 5 6
0
2
4
6
8
10
12
x
x (2+sin(2 x))
Adaptive Quadrature
? ?10 2 31dxx
4
1
2
11 2 ??
?
??
?
??M
2
1
2
101 2 ???T
12
1
4
1
3
1 ??
6
1
2
1
3
1 ???
)(2)( MSTS ????
TMS 3132 ??
Simpson’s rule:
2) ),()(4)((6
bacbfcfafhS ?????
TMS 3132 ??
四阶求积公式
Composite quadrature rule,对 [a,c],[c,d],d=(a-c)/2,e=(c-b)/2
)),()(4)(2)(4)((122 bfefcfdfafhS ?????
Adaptive Quadrature
S和 S2都逼近同样的积分，他们的差可以用来估计误差
)( 2 SSE ??
二者也可以组合起来得到更精确的逼近。这两个求积
公式都是四阶的，但 S2的步长是 S的一半
)(16)( 2SQSQ ???
15/)( 22 SSSQ ???
六阶求积公式
Weddle’s rule,the sixth order Newton-Cotes rule or the
first step of Romberg integration.
Extrapolated Simpson’s rule
quadtx
Matlab中 quad用递归自适应的外插 Simpon公式
quadtx是一个简化版本。先求出三点的函数值，然后用递
归过程 quadtxstep完成积分。
quadtx.m
[Qa,ka] = quadtxstep(F,a,c,tol,fa,fd,fc,varargin{:});
[Qa,ka] = quadtxstep(F,c,b,tol,fc,fe,fb,varargin{:});
tol的故事：
http:// www.inf.ethz.ch/personal/gander,
Specifying Integrands
? ?10 41 1 dxx
f = inline(‘1/sqrt(1+x^4)’)
Q = quadtx(f,0,1)
?10 )s i n ( dxx x
f = inline(‘sin(x)/x’)
Q = quadtx(f,0,pi)
Q = quadtx(f,realmin,pi)
Q = quadtx(@sinc,0,pi)
? ?? ?10 11 )1( dxtt z ? f = inline(‘t^(z-1)*(1-t)^(w-1)’,’t’,’z’,’w’)
z = 8/3; w = 10/3; tol = 1.e-6;
Q = quadtx(f,0,1,tol,z,w)
vectorize
10
20
30
40
50
60
70
80
90
100
93
Performance
04.0)9.0(
1
01.0)3.0(
1)(
22 ?????? xxxh
humps 函数
quadgui(@humps,0,1,1.e-4)
Performance
syms x
h = 1/((x-.3)^2+.01) + 1/((x-.9)^2+.04) - 6
I = int(h)
D = simple(int(h,0,1))
Qexact = double(D)
for k=1:12
tol=10^(-k); [Q,fcount]=quadtx(@humps,0,1,tol);
err=Q-Qexact; ratio=err/tol;
fprintf('%8.0e %21.14f %7d %13.3e %9.3f\n',...
tol,Q,fcount,err,ratio),end
tol Q count err ratio
1e-001 29.83328444174864 25 -2.504e-002 -0.250
1e-002 29.85791444629948 41 -4.109e-004 -0.041
1e-003 29.85834299237637 69 1.760e-005 0.018
1e-004 29.85832444437543 93 -9.511e-007 -0.010
1e-005 29.85832551548643 149 1.200e-007 0.012
1e-006 29.85832540194041 265 6.442e-009 0.006
1e-007 29.85832539499819 369 -5.005e-010 -0.005
1e-008 29.85832539552631 605 2.763e-011 0.003
1e-009 29.85832539549604 1061 -2.640e-012 -0.003
1e-010 29.85832539549890 1469 2.274e-013 0.002
1e-011 29.85832539549867 2429 -7.105e-015 -0.001
1e-012 29.85832539549867 4245 0.000e+000 0.000
