x0 x1 x2 x3 x
lim y (k ) n 1
y n1
k
11
向后欧拉公式的误差
yn1 yn hf (xn1, yn1), n 0, 1,
Tn1
y(xn1)
yn1
h2 2
y(xn) O(h3)
向前欧拉公式 局部误差主项
h2 2
y ( x n )
~QA
向后欧拉公式 局部误差主项
h2 2
y ( x n )
~QB
y
向前、向后欧拉公式的右端平
均，则两个误差主项刚好抵消
Pn
向后欧拉 公式具有 1阶精度
A Q B
梯形公式
xn xn+1
x
12
向前欧拉公式
yn1 yn hf (xn, yn)
向后欧拉公式
yn1 yn hf (xn1, yn1)
梯形 公式
yn1
yn
h 2
[
f
( xn ,
di i (1 i )
dt i ( 0 ) i 0
Logistic 模型
i(t)
1
1
1 i0
1 e t
3
传染病模型
传染病有免疫性——病人治愈 后即移出感染系统，称移出者
SIR模型
假设 1）总人数N不变，病人、健康人和移
出者的比例分别为 i ( t ), s ( t ), r ( t )
Tn 1
y( xn1)
y n 1
h2 2
y(xn ) O(h3 ) O(h2 )
若一种算法的局部截断误差为 O(h p1 ) ，则称该算法
具有 p 阶精度
向前欧拉公式具有1阶精度
局部截断误差主项为
h2 2
y ( x n )
10
欧拉 方法
x n 1
y ( xn 1 ) y ( xn ) f ( x , y ( x )) d x,
梯形公式具有2阶精度 局部截断 误差主项
h3
12
y'''(xn )
隐式公式迭代求解
y(0) n1
yn
hf
(xn ,
yn )
y (k 1) n 1
yn
h[f 2
(xn , yn )
f
(
xn 1 ,
y
(k ) n 1
)]
k 0,1,2,, n 0,1,2,
xn1
y(xn1) y(xn) f (x, y(x))dx
xn
取不同的点
各种欧拉公式
y y0 P0
P1 P2 P3 y=y(x)
x0 x1 x2 x3 x
8
向前欧拉公式
x取左端点xn
y ( x n 1 ) y ( x n ) hf ( x n , y ( x n ))
近似：yny(xn) ,yn+1 y(xn+1)
常微分方程数值解
主要内容
1. 两个最常用的数值解法: • 欧拉（Euler）方法 • 龙格-库塔（Runge-Kutta）方法
2. 龙格-库塔方法的MATLAB实现
1
传染病模型
SI 模型 区分已感染者(病人)和未感染者(健康人)
假设 1）总人数N不变，病人和健康
人的 比例分别为 i(t ), s(t )
2）每个病人每天有效接触人数 ~ 日
为, 且使接触的健康人致病
接触率
建模 N[i(t t) i(t)] [s(t)]Ni(t)t
di si
dt
s(t) i(t) 1
di i (1 i )
dt i ( 0 ) i 0
2
传染病模型
SI 模型 区分已感染者(病人)和未感染者(健康人)
以推广到解常 微分方程组、 高阶微分方程
2阶龙格—库塔公式 具有2阶精度
14
龙格-库塔方法
y f ( x, y ), y ( x0 ) y0
xn
向后欧拉公式 x取右端点xn+1,
yn1 yn hf (xn1, yn1), n 0, 1, 右端yn+1未知
隐式公式，迭代求解
y P0
y0
P2 P1
y=y(x) P3
y(0) n1
yn
hf
(xn ,
yn )
y ( k 1) n1
yn
hf
( xn1 ,
y(k) n1
)
k 0,1,2,, n 0,1,2,
di
dt
si i
ds
dt
si
i(0)
i0 ,
s(0)
s0
无法求出 i(t ), s (t )
的解析解
i0 s0 1 (通常r(0) r0很小）
5
相平面分析方法（Phase Plane Analysis)
di
dt
si i
s(t),i(t) 时间t的相（像）
ds
2）病人的日接触率 , 日治愈率, 接触数 = /
建模 s (t ) i(t ) r (t ) 1
需建立 i ( t ), s ( t ), r ( t ) 的两个方程
4
传染病模型
SIR模型
N[i(t t) i(t)] Ns(t)i(t)t Ni(t)t
N[s(t t) s(t)] Ns(t)i(t)t
13
改进的欧拉公式
将梯形公式的迭代过程简化为两步
yn1 yn hf (xn , yn )
预测
yn 1
yn
h[ 2
f
( xn ,
yn)
f
( xn 1,
yn1)]
校正
n 0,1,2,
可
y n 1
yn
h 2
(k1
k2)
欧拉公式都可
表 示
k1 f (xn , yn )
为
k 2 f ( x n1 , y n hk 1 )
yn1 yn hf (xn, yn), n 0,1,
9
误差 分析
假设到第n步公式右端yn没有误差,即yn=y(xn), 从xn到xn+1一步的计算值yn+1与精确值y(xn+1) 之差,称为局部截断误差.
y ( xn1 )
y(xn ) hy(xn )
h2 2
y(xn ) O(h3)
y n1 y ( xn ) hf ( xn , y ( xn )) y ( xn ) hy ( xn )
不求解析解 y y(x) ，而在一系列离散点
x0 x1 x2 xn
求y(xn )的近似值 y
记作yn (n 1,2,)
y y(x)
通常取等步长h
xn x0 nh
y0
y1 y2
x0 x1 x 2
yn
xn x
7
欧拉方法
y f ( x, y), y( x0 ) y0
基本思路
si
dt
s i平面称为相平面
i(0)
i0 ,
s(0)
s0
(s(t), i(t)) 的轨迹称为相轨线
相平面分析是在相平面上分析相轨线随时间变 化的变化情况，是一种定性分析的方法。它也 是应用中非常重要的一种分析方法。
6
“常微分方程初值问题数值解”的提法
设 y f ( x, y ), y ( x0 ) y0的解 y y ( x)存在且唯一