Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
常微分方程
Integrating Differential Equatioins
00 )()),(,(
)( ytytytf
dt
tdy ??
,.,,1,0),( ?? ntyy nn
nnn tth ?? ?1
? ???? htt dssysftyhty ))(,()()(
? ???? htt dssftyhty )()()(
特殊情况：
积分方程：
步长,自动确定
微分方程：
数值解：
System of Equations
)(txx ????
??
?
??
??
)(
)()(
tx
txty
? ??
?
??
?
????
?
??
?
?? )(
)(
)(
)()(
1
2
ty
ty
tx
txty ??
function ydot = harmonic(t,y)
ydot=[y(2); -y(1)]
y=inline(‘[0 1;-1 0]*y’,’t’,’y’);
Two Body Problem
3
3
)(/)()(
)(/)()(
trtvtv
trtutu
??
??
??
?? 22 )()()( tvtutr ??
?
?
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
)(
tv
tu
tv
tu
ty
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
3
)(/)(
)(/)(
)(
)(
)(
trtv
trtu
tv
tu
ty
?
?
?
function ydot =twobody(t,y)
r=sqrt(y(1)^2+y(2)^2);
ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3];
%ydot=[y(3:4;-y(1:2)/norm(y(1:2))^3];
Linearized Differential Equations
)()(),(),( cccc yyJttytfytf ????? ?
),(),,( cccc ytyfJyttf ??????? ??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
),.,,,,(
),.,,,,(
),.,,,,(
)(
)(
)(
1
12
11
2
1
nn
n
n
n
yytf
yytf
yytf
ty
ty
ty
dt
d
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
n
nnn
n
n
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
J
?
????
?
?
21
2
2
2
1
2
1
2
1
1
1
Jyy ??
Linearized Differential Equations
J的特征值是
kkk i??? ?? )(d iag k???
1??? VVJ yVx ?
kkk xx ??? )()( )( cttk txetx ck ?? ?
0?k?
0?k?
0?k?
解增长
解衰减
解振荡
Linearized Differential Equations
yy ?
?
?
??
?
?? 01
10?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
2
3
1
4
3
)(/)(
)(/)(
)(
)(
trty
trty
ty
ty
y?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0023
0032
000
000
1
2
1
2
221
21
2
2
2
1
5
5
5
yyyy
yyyy
r
r
r
J
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
i
i
r 2
2
1
2/3
?
例 1：
例 2：
itit ee ?,
Single Step Methods
Euler 方法
htt
ythfyy
nn
nnnn
??
??
?
?
1
1 ),(
t=t0; y=y0;
while t<=tfinal
y = y+h*feval(f,t,y)
t=t+h
end
误差,O(h),需要小步长；不能自动确定步长
Single Step Methods
htt
hsyy
s
h
y
h
tfs
ytfs
nn
nn
nn
nn
??
??
???
?
?
?
1
21
12
1
)
2
,
2
(
),(
htt
ss
hyy
hsyhtfs
ytfs
nn
nn
nn
nn
??
?
??
???
?
?
?
1
21
1
12
1
2
),(
),(midpoint trapezoid
如果同时用两个方法计算，则在一点得到不同的值，
这样 我们可以得到一个误差估计，并用来选取步长。
另外，这两个值的外插组合可以得到更精确的结果。
Single Step Methods
单步方法：积分思想。
? ???? htt dssysftyhty ))(,()()(
,),())(,( ?? ?? ? i
i
iii
ht
t
yysfchdssysf
算法关键：自动确定时间步长 h
htt
ssss
h
yy
hsyhtfs
s
h
y
h
tfs
s
h
y
h
tfs
ytfs
nn
nn
nn
nn
nn
nn
??
?????
???
???
???
?
?
?
1
43211
34
23
12
1
)22(
6
),(
)
2
,
2
(
)
2
,
2
(
),(
Single Step Methods
Runge-Kutta方法,1905（德国）
若 f不依赖 y,则为
simpson公式
如何确定 h?
Single Step Methods
??
?
????
1
1
,,.,,,1),,(
i
j
jjinini kishyhtfs ??
?
?
? ??
k
i
iinn shyy
1
1 ?
?
?
? ?
k
i
iin she
1
1 ?
两组方法
计算方法：
如果误差小于指定的精度，进行下一步，否则，缩小步
长。
Matlab,ode23,ode45
BS23 algorithm
ode23,Bogacki,Shampine(1989)和 Shampine(1994),”23”表
示用两算法：一个 2阶，一个 3阶
)9865(
72
),(
)432(
9
)
4
3
,
4
3
(
)
2
,
2
(
),(
43211
114
3211
1
23
12
1
ssss
h
e
ytfs
sss
h
yy
htt
hsyhtfs
s
h
y
h
tfs
ytfs
n
nn
nn
nn
nn
nn
nn
?????
?
????
??
???
???
?
?
??
?
?
ode23tx.m
Bogacki,P,and LF
Shampine,"A 3(2) pair
of Runge-Kutta
formulas," Appl.
Math,Letters,Vol,2,
1989,pp 1-9.
Examples
ode23tx(‘0’,[0 10],1)
ode23tx(‘2*y-y^2’,[0 10],1) 2/(1+exp(-2*t))
1
精确解
F=inline('[y(2);-y(1)]','t','y')
ode23tx(F,[0 2*pi],[1;0])
[t,y]=ode23tx(F,[0 2*pi],[1;0])
plot(y(:,1),y(:,2),’-o’)
axis([-1.2 1.2 –1.2 1.2])
axis square
opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’,@odephas2)
Examples
ode23tx(@twobody,[0 2*pi],[1;0;0;1]);
0 1 2 3 4 5 6 7-1.5
-1
-0.5
0
0.5
1
1.5
Examples
y0=[1;0;0;3]; ode23(@twobody,[0 2*pi],y0);
0 1 2 3 4 5 6 7-2
0
2
4
6
8
10
12
14
16
18
Examples
y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0);
plot(y(:,1),y(:,2),’-’,0,0,’ro’); axis equal
- 1, 5 -1 - 0, 5 0 0, 5 1
0
2
4
6
8
10
12
14
16
18
Lorenz Attractor
Ayy ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
1
0
0
2
2
?
??
?
y
y
A
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
3
2
1
ty
ty
ty
ty
y1大气对流,y2,y3为水平和垂直的温度变化,sigma是
Prandtl数,rho是规范 Rayleigh数,beta依赖区域。
3/8,28,10 ??? ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
1
0
0
??
??
??
A
A是奇异的，当且仅当
)1( ??? ???
?
?
?
?
?
?
?
?
?
? ?
?
?
?
? 1
0y
00 )( yty ?
?
?
?
?
?
?
?
?
?
?
?
0
0
0
)(ty?
这两个点是不稳定的
天气预报的准确性：
http://www.nani.com.tw/senior/seniorteach/stfastnews/stnfast
news/stnfastnews4/stnfastnews4_2.htm
Lorenz现象的数学：
http://www.sciscape.org/news_detail.php?news_id=225
分形艺术电子版：
http://www.phil.pku.edu.cn/personal/huajie/fractalart/html/no
tice.htm
混沌理论：
http://www.lamost.org/~yzhao/doc/chaos.html
rho = 28; sigma = 10; beta = 8/3;
eta = sqrt(beta*(rho-1));
A = [ -beta 0 eta
0 -sigma sigma
-eta rho -1 ];
yc = [rho-1; eta; eta];
y0 = yc + [0; 0; 3];
tspan = [0 Inf];
opts = odeset('reltol',3.e-4,'outputfcn',@lorenzplot);
ode23(@lorenzeqn,tspan,y0,opts,A);
气象学家 Lorenz提出一篇论文，名叫「一只蝴蝶拍一下翅膀会不会
在 Taxas州引起龙卷风？」论述某系统如果初期条件差一点点，结果
会很不稳定，他把这种现象戏称做「 蝴蝶效应 」。就像我们投掷骰
子两次，无论我们如何刻意去投掷，两次的物理现象和投出的点数
也不一定是相同的。 Lorenz为何要写这篇论文呢？
这故事发生在 1961年的某个冬天，他如往常一般在办公室操作气象
电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑
就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，
因此模拟出气象变化图。
这一天,Lorenz想更进一步了解某段纪录的後续变化，他把某时刻
的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，
电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡
并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆
。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差
异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题
是他输入的数据差了 0.000127，而这些微的差异却造成天壤之别。
所以长期的准确预测天气是不可能的。
0 5 10 15 20 25 30
- 3 0
- 2 0
- 1 0
0
10
20
30
40
50
0 5 10 15 20 25 30-80
-60
-40
-20
0
20
40
60
y1
y2
y3
5
10 15
20
25 30
35 40
45
-30
-20
-10
0
10
-80
-70
-60
-50
-40
-30
-20
plot3(y(:,1),y(:,2),y(:,3))
0 5 10 15 20 25 30-30
-20
-10
0
10
20
30
40
50
3e.-10
Stiffness
A problem is stiff if the solution being sought is varying
slowly,but there are nearby solutions that very rapidly,so the
numerical method must take small steps to obtain satisfactory
results.
A model of flame propagation:
?
?
/20
)0(
32
??
?
??
t
y
yyy? y是球的半径,y^2和
y^3与球的表面积和体
积有关
想一下，点燃一根火柴，光球迅速增长，到达一个关键
的大小，然后维持它的大小（由于进入球内氧气和消耗
的氧气平衡）
0 10 20 30 40 50 60 70 80 90 1000
0.2
0.4
0.6
0.8
1
1.2
1.4
eta=0.02;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
5.03 5.035 5.04 5.045 5.05 5.055 5.06
x 104
0.9985
0.999
0.9995
1
1.0005
1.001
1.0015
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
5, 0 2 7 5 5, 0 2 8 5, 0 2 8 5 5, 0 2 9 5, 0 2 9 5 5, 0 3 5, 0 3 0 5 5, 0 3 1 5, 0 3 1 5 5, 0 3 2 5, 0 3 2 5
x 1 0
4
0, 9 9 8
0, 9 9 8 5
0, 9 9 9
0, 9 9 9 5
1
1, 0 0 0 5
1, 0 0 1
t o l = 1, e - 6
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta,opts);
0 1 2 3 4 5 6 7 8 9 10
x 104
0
0.2
0.4
0.6
0.8
1
1.2
1.4
ode23s
eta=0.00002; ode23s(inline('y^2-y^3','t','y'),[0 2/eta],eta);
Events
如何确定计算的终结条件，是一件重要的事情。例如物
体下落问题中物体何时落地？两体问题中，轨道的周期？
0))(,(
)(
),(
**
00
?
?
?
tytg
yty
ytfy? 需要有一个终结条件 g
物体下落问题
0)(
0)0(,1)0(
1
*
2
?
??
???
ty
yy
yy
?
???
function ydot=f(t,y)
ydot=[y(2); -1+y(2)^2];
Events
function [gstop,isterminal,directioin]=g(t,y)
gstop=y(1); isterminal=1; direction=[];
opts=odeset('events',@g)
y0=[1;0];
[t,y,tfinal]=ode45(@f,[0 Inf],y0,opts); tfinal
plot(t,y(:,1),'-',[0 tfinal],[1 0],'o')
axis([-.1 tfinal+.1 -.1 1.1])
xlabel('t'); ylabel('y');title('Falling body')
text(1.2,0,['tfinal=' num2str(tfinal)])
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0
0.2
0.4
0.6
0.8
1
t
y
Falling body
tfinal=1.6575
自由落体
0 0.2 0.4 0.6 0.8 1
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
两体问题
orbit.m
twobody.m
gstop.m
多步方法
单步方法,内存，少
后一点的值，由前面一点导出
合适的步长选取
多步方法,内存，多
后一点的值，由前面多点导出
合适的步长选取
对于光滑的解，多步方法比单步方法精确
Matlab's ODE Solvers
非刚性系统：
ode45(Runge-Kutta45)
ode23(Runge-Kuatta23)
ode113(Adams-Bashforth-Moulton PECE),多步方法
刚性系统：
ode15s(Gear方法 ),多步方法
ode23s(二阶 modified Rosenbroack formula),单步
ode23t(trapezoidal rule),solve DAEs
ode23tb(TR-BDF2)
Errors
舍入误差：浮点表示
离散误差：局部离散误差，全局离散误差
? ?
?
?Nt
t
N
nn tfhdf
0
1
0
)()( ??
局部离散误差
? ??? 1 )()( nnttnnn dftfhd ??
全局离散误差 ?? ?? ? Nt
t
N
nnN dftfhe
0
)()(
1
0
??
?? tt dfty 0 )()( ??
一般的方程,f(t,y)在局部离散误差假设前一时刻解是精
确的，估计下一时刻的误差。
Errors
1|| ?? pnn Chd一个方法是 p阶 (Order):
Euler方法,),(
1 nnnnn ytfhyy ???
)()( 211 nnnnn hOtuyd ??? ??局部离散误差
全局离散误差 )()(* 2 hOhONe
n ??
步长缩小一半，局部误差,1/4，单需要两步达到原来的
一步，全局误差,1/2.
可以用对 f(t,y)是一个不依赖 y的多项式来检测方法
的阶：如果对 t的 p-1次是精确的，而对 t的 p次是不精
确的，则方法是 p阶的。
Errors
假设每一步的舍入误差是 epsilon,N(L/h)步后的影响是
h
LN ?
? ?
考虑离散误差的影响，总误差是
h
LCh p ??
如果两者的误差是可比的，则
1
1
?
?
?
??
?
?? p
C
Lh ?
1
1
?
?
?
??
?
?? p
L
CLN
?
p N
1 4.5*10^17
3 5647721
5 37285
10 864
Euler 方法
ode23
ode45
ode113
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
常微分方程
Integrating Differential Equatioins
00 )()),(,(
)( ytytytf
dt
tdy ??
,.,,1,0),( ?? ntyy nn
nnn tth ?? ?1
? ???? htt dssysftyhty ))(,()()(
? ???? htt dssftyhty )()()(
特殊情况：
积分方程：
步长,自动确定
微分方程：
数值解：
System of Equations
)(txx ????
??
?
??
??
)(
)()(
tx
txty
? ??
?
??
?
????
?
??
?
?? )(
)(
)(
)()(
1
2
ty
ty
tx
txty ??
function ydot = harmonic(t,y)
ydot=[y(2); -y(1)]
y=inline(‘[0 1;-1 0]*y’,’t’,’y’);
Two Body Problem
3
3
)(/)()(
)(/)()(
trtvtv
trtutu
??
??
??
?? 22 )()()( tvtutr ??
?
?
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
)(
tv
tu
tv
tu
ty
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
3
)(/)(
)(/)(
)(
)(
)(
trtv
trtu
tv
tu
ty
?
?
?
function ydot =twobody(t,y)
r=sqrt(y(1)^2+y(2)^2);
ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3];
%ydot=[y(3:4;-y(1:2)/norm(y(1:2))^3];
Linearized Differential Equations
)()(),(),( cccc yyJttytfytf ????? ?
),(),,( cccc ytyfJyttf ??????? ??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
),.,,,,(
),.,,,,(
),.,,,,(
)(
)(
)(
1
12
11
2
1
nn
n
n
n
yytf
yytf
yytf
ty
ty
ty
dt
d
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
n
nnn
n
n
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
J
?
????
?
?
21
2
2
2
1
2
1
2
1
1
1
Jyy ??
Linearized Differential Equations
J的特征值是
kkk i??? ?? )(d iag k???
1??? VVJ yVx ?
kkk xx ??? )()( )( cttk txetx ck ?? ?
0?k?
0?k?
0?k?
解增长
解衰减
解振荡
Linearized Differential Equations
yy ?
?
?
??
?
?? 01
10?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
2
3
1
4
3
)(/)(
)(/)(
)(
)(
trty
trty
ty
ty
y?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0023
0032
000
000
1
2
1
2
221
21
2
2
2
1
5
5
5
yyyy
yyyy
r
r
r
J
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
i
i
r 2
2
1
2/3
?
例 1：
例 2：
itit ee ?,
Single Step Methods
Euler 方法
htt
ythfyy
nn
nnnn
??
??
?
?
1
1 ),(
t=t0; y=y0;
while t<=tfinal
y = y+h*feval(f,t,y)
t=t+h
end
误差,O(h),需要小步长；不能自动确定步长
Single Step Methods
htt
hsyy
s
h
y
h
tfs
ytfs
nn
nn
nn
nn
??
??
???
?
?
?
1
21
12
1
)
2
,
2
(
),(
htt
ss
hyy
hsyhtfs
ytfs
nn
nn
nn
nn
??
?
??
???
?
?
?
1
21
1
12
1
2
),(
),(midpoint trapezoid
如果同时用两个方法计算，则在一点得到不同的值，
这样 我们可以得到一个误差估计，并用来选取步长。
另外，这两个值的外插组合可以得到更精确的结果。
Single Step Methods
单步方法：积分思想。
? ???? htt dssysftyhty ))(,()()(
,),())(,( ?? ?? ? i
i
iii
ht
t
yysfchdssysf
算法关键：自动确定时间步长 h
htt
ssss
h
yy
hsyhtfs
s
h
y
h
tfs
s
h
y
h
tfs
ytfs
nn
nn
nn
nn
nn
nn
??
?????
???
???
???
?
?
?
1
43211
34
23
12
1
)22(
6
),(
)
2
,
2
(
)
2
,
2
(
),(
Single Step Methods
Runge-Kutta方法,1905（德国）
若 f不依赖 y,则为
simpson公式
如何确定 h?
Single Step Methods
??
?
????
1
1
,,.,,,1),,(
i
j
jjinini kishyhtfs ??
?
?
? ??
k
i
iinn shyy
1
1 ?
?
?
? ?
k
i
iin she
1
1 ?
两组方法
计算方法：
如果误差小于指定的精度，进行下一步，否则，缩小步
长。
Matlab,ode23,ode45
BS23 algorithm
ode23,Bogacki,Shampine(1989)和 Shampine(1994),”23”表
示用两算法：一个 2阶，一个 3阶
)9865(
72
),(
)432(
9
)
4
3
,
4
3
(
)
2
,
2
(
),(
43211
114
3211
1
23
12
1
ssss
h
e
ytfs
sss
h
yy
htt
hsyhtfs
s
h
y
h
tfs
ytfs
n
nn
nn
nn
nn
nn
nn
?????
?
????
??
???
???
?
?
??
?
?
ode23tx.m
Bogacki,P,and LF
Shampine,"A 3(2) pair
of Runge-Kutta
formulas," Appl.
Math,Letters,Vol,2,
1989,pp 1-9.
Examples
ode23tx(‘0’,[0 10],1)
ode23tx(‘2*y-y^2’,[0 10],1) 2/(1+exp(-2*t))
1
精确解
F=inline('[y(2);-y(1)]','t','y')
ode23tx(F,[0 2*pi],[1;0])
[t,y]=ode23tx(F,[0 2*pi],[1;0])
plot(y(:,1),y(:,2),’-o’)
axis([-1.2 1.2 –1.2 1.2])
axis square
opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’,@odephas2)
Examples
ode23tx(@twobody,[0 2*pi],[1;0;0;1]);
0 1 2 3 4 5 6 7-1.5
-1
-0.5
0
0.5
1
1.5
Examples
y0=[1;0;0;3]; ode23(@twobody,[0 2*pi],y0);
0 1 2 3 4 5 6 7-2
0
2
4
6
8
10
12
14
16
18
Examples
y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0);
plot(y(:,1),y(:,2),’-’,0,0,’ro’); axis equal
- 1, 5 -1 - 0, 5 0 0, 5 1
0
2
4
6
8
10
12
14
16
18
Lorenz Attractor
Ayy ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
1
0
0
2
2
?
??
?
y
y
A
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
3
2
1
ty
ty
ty
ty
y1大气对流,y2,y3为水平和垂直的温度变化,sigma是
Prandtl数,rho是规范 Rayleigh数,beta依赖区域。
3/8,28,10 ??? ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
1
0
0
??
??
??
A
A是奇异的，当且仅当
)1( ??? ???
?
?
?
?
?
?
?
?
?
? ?
?
?
?
? 1
0y
00 )( yty ?
?
?
?
?
?
?
?
?
?
?
?
0
0
0
)(ty?
这两个点是不稳定的
天气预报的准确性：
http://www.nani.com.tw/senior/seniorteach/stfastnews/stnfast
news/stnfastnews4/stnfastnews4_2.htm
Lorenz现象的数学：
http://www.sciscape.org/news_detail.php?news_id=225
分形艺术电子版：
http://www.phil.pku.edu.cn/personal/huajie/fractalart/html/no
tice.htm
混沌理论：
http://www.lamost.org/~yzhao/doc/chaos.html
rho = 28; sigma = 10; beta = 8/3;
eta = sqrt(beta*(rho-1));
A = [ -beta 0 eta
0 -sigma sigma
-eta rho -1 ];
yc = [rho-1; eta; eta];
y0 = yc + [0; 0; 3];
tspan = [0 Inf];
opts = odeset('reltol',3.e-4,'outputfcn',@lorenzplot);
ode23(@lorenzeqn,tspan,y0,opts,A);
气象学家 Lorenz提出一篇论文，名叫「一只蝴蝶拍一下翅膀会不会
在 Taxas州引起龙卷风？」论述某系统如果初期条件差一点点，结果
会很不稳定，他把这种现象戏称做「 蝴蝶效应 」。就像我们投掷骰
子两次，无论我们如何刻意去投掷，两次的物理现象和投出的点数
也不一定是相同的。 Lorenz为何要写这篇论文呢？
这故事发生在 1961年的某个冬天，他如往常一般在办公室操作气象
电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑
就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，
因此模拟出气象变化图。
这一天,Lorenz想更进一步了解某段纪录的後续变化，他把某时刻
的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，
电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡
并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆
。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差
异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题
是他输入的数据差了 0.000127，而这些微的差异却造成天壤之别。
所以长期的准确预测天气是不可能的。
0 5 10 15 20 25 30
- 3 0
- 2 0
- 1 0
0
10
20
30
40
50
0 5 10 15 20 25 30-80
-60
-40
-20
0
20
40
60
y1
y2
y3
5
10 15
20
25 30
35 40
45
-30
-20
-10
0
10
-80
-70
-60
-50
-40
-30
-20
plot3(y(:,1),y(:,2),y(:,3))
0 5 10 15 20 25 30-30
-20
-10
0
10
20
30
40
50
3e.-10
Stiffness
A problem is stiff if the solution being sought is varying
slowly,but there are nearby solutions that very rapidly,so the
numerical method must take small steps to obtain satisfactory
results.
A model of flame propagation:
?
?
/20
)0(
32
??
?
??
t
y
yyy? y是球的半径,y^2和
y^3与球的表面积和体
积有关
想一下，点燃一根火柴，光球迅速增长，到达一个关键
的大小，然后维持它的大小（由于进入球内氧气和消耗
的氧气平衡）
0 10 20 30 40 50 60 70 80 90 1000
0.2
0.4
0.6
0.8
1
1.2
1.4
eta=0.02;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
5.03 5.035 5.04 5.045 5.05 5.055 5.06
x 104
0.9985
0.999
0.9995
1
1.0005
1.001
1.0015
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta);
5, 0 2 7 5 5, 0 2 8 5, 0 2 8 5 5, 0 2 9 5, 0 2 9 5 5, 0 3 5, 0 3 0 5 5, 0 3 1 5, 0 3 1 5 5, 0 3 2 5, 0 3 2 5
x 1 0
4
0, 9 9 8
0, 9 9 8 5
0, 9 9 9
0, 9 9 9 5
1
1, 0 0 0 5
1, 0 0 1
t o l = 1, e - 6
eta=0.00002;ode23tx(‘y^2-y^3’,[0 2/eta],eta,opts);
0 1 2 3 4 5 6 7 8 9 10
x 104
0
0.2
0.4
0.6
0.8
1
1.2
1.4
ode23s
eta=0.00002; ode23s(inline('y^2-y^3','t','y'),[0 2/eta],eta);
Events
如何确定计算的终结条件，是一件重要的事情。例如物
体下落问题中物体何时落地？两体问题中，轨道的周期？
0))(,(
)(
),(
**
00
?
?
?
tytg
yty
ytfy? 需要有一个终结条件 g
物体下落问题
0)(
0)0(,1)0(
1
*
2
?
??
???
ty
yy
yy
?
???
function ydot=f(t,y)
ydot=[y(2); -1+y(2)^2];
Events
function [gstop,isterminal,directioin]=g(t,y)
gstop=y(1); isterminal=1; direction=[];
opts=odeset('events',@g)
y0=[1;0];
[t,y,tfinal]=ode45(@f,[0 Inf],y0,opts); tfinal
plot(t,y(:,1),'-',[0 tfinal],[1 0],'o')
axis([-.1 tfinal+.1 -.1 1.1])
xlabel('t'); ylabel('y');title('Falling body')
text(1.2,0,['tfinal=' num2str(tfinal)])
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0
0.2
0.4
0.6
0.8
1
t
y
Falling body
tfinal=1.6575
自由落体
0 0.2 0.4 0.6 0.8 1
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
两体问题
orbit.m
twobody.m
gstop.m
多步方法
单步方法,内存，少
后一点的值，由前面一点导出
合适的步长选取
多步方法,内存，多
后一点的值，由前面多点导出
合适的步长选取
对于光滑的解，多步方法比单步方法精确
Matlab's ODE Solvers
非刚性系统：
ode45(Runge-Kutta45)
ode23(Runge-Kuatta23)
ode113(Adams-Bashforth-Moulton PECE),多步方法
刚性系统：
ode15s(Gear方法 ),多步方法
ode23s(二阶 modified Rosenbroack formula),单步
ode23t(trapezoidal rule),solve DAEs
ode23tb(TR-BDF2)
Errors
舍入误差：浮点表示
离散误差：局部离散误差，全局离散误差
? ?
?
?Nt
t
N
nn tfhdf
0
1
0
)()( ??
局部离散误差
? ??? 1 )()( nnttnnn dftfhd ??
全局离散误差 ?? ?? ? Nt
t
N
nnN dftfhe
0
)()(
1
0
??
?? tt dfty 0 )()( ??
一般的方程,f(t,y)在局部离散误差假设前一时刻解是精
确的，估计下一时刻的误差。
Errors
1|| ?? pnn Chd一个方法是 p阶 (Order):
Euler方法,),(
1 nnnnn ytfhyy ???
)()( 211 nnnnn hOtuyd ??? ??局部离散误差
全局离散误差 )()(* 2 hOhONe
n ??
步长缩小一半，局部误差,1/4，单需要两步达到原来的
一步，全局误差,1/2.
可以用对 f(t,y)是一个不依赖 y的多项式来检测方法
的阶：如果对 t的 p-1次是精确的，而对 t的 p次是不精
确的，则方法是 p阶的。
Errors
假设每一步的舍入误差是 epsilon,N(L/h)步后的影响是
h
LN ?
? ?
考虑离散误差的影响，总误差是
h
LCh p ??
如果两者的误差是可比的，则
1
1
?
?
?
??
?
?? p
C
Lh ?
1
1
?
?
?
??
?
?? p
L
CLN
?
p N
1 4.5*10^17
3 5647721
5 37285
10 864
Euler 方法
ode23
ode45
ode113
