5
3、使用泰勒公式
以此方法为基础，有龙格-库塔法、线性多步法等方法
库塔三阶方法 四阶龙格-库塔公式
? ?
yn ? 1
?
?
yn ?
h 6
(K1
?
4K2 ?
K 3 ),
? ?
K
1
?
f ( x n , y n ),
? ?
K
2
?
f ( xn
?
h 2 , yn
?
h 2 K 1 ),
?? K 3 ? f ( x n ? h , y n ? hK 1 ? 2 hK 2 ).
6
4、数值公式的精度 当一个数值公式的截断误差可表示为o(hk)时(k为正整数，h为步长)，称它是一个k
阶公式。 k越大，则数值公式的精度越高。
欧拉法是一阶公式，改进的欧拉法是二阶公式。
线性多步法有四阶阿达姆斯外插公式和内插公式。
北京邮电大学数学系
7
(三)用Matlab 软件求常微分方程的数值解 [t，x]= solver( ‘f', ts, x0, options ）
[t,y]=ode45(@dfun2,[0,20],[2,0]);
1、建立m-文件dfun2.m 如下： function dx=dfun2(t,y)
)]
k ? 0,1,2,?
对于已给的精确度
?，当满足
y ? y ( k ? 1)
(k)
i?1
i?1
?
? 时，取
y i?1
?
y ， （ k ?1） i?1
然后继续下一步 y i? 2 的计算。
此即改进的欧拉法
北京邮电大学数学系
4
? 中点欧拉公式 /* midpoint formula */
中心差商近似导数
北京邮电大学数学系
2
(二)建立数值解法的一些途径
设 x i?1 ? xi ? h, i ? 0,1,2,? n ? 1, 可用以下离散化方法求
? y' ? f(x, y)
? ?
y(x
0
)
?
y0
1、用差商代替导数
解微分方程：
若步长h较小，则有
故有公式：
y'( x) ? y( x ? h) ? y( x) h
第五章 微分方程模型
常微分方程的数值解 5.1 传染病模型
北京邮电大学数学系
1
常微分方程的数值解及实验
(一)常微分方程数值解
高数中微分方程解法在实际中基本不会直接使用
在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题， 一般是要求得到解在若干个点上满足规定精确度的近似值 ，或者得到一个满足精确度要求的便 于计算的表达式。
y?( x1 ) ?
y( x2 ) ? y( x0 ) 2h
y￠( x1 )
y( x 2 ) ? y( x 0 ) ? 2 h f ( x 1 , y( x 1 ))
x0
x1
x2
yi ? 1 ? yi ? 1 ? 2 h f ( x i , yi ) i ? 1, ... , n ? 1
北京邮电大学数学系
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令 为：options=odeset('reltol',rt,'abstol',at), rt, at：分别为设定的相对误差和绝对误差.
北京邮电大学数学系
8
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分 量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.
北京邮电大学数学系
9
微分方程求解实例
设取步长
，从
到
用四阶龙格-库塔方法
求解初值问题
建立m-文件dfun1.m如下
function dx=dfun1(x,y) dx=y-2*x/y;
输入命令
h=0.2; ts=0:h:1; y0=1;
[t, x]=ode45('dfun1',ts,y0);
[t, x], plot(t,x)
北京邮电大学数学系
3、结果如图
0
1.0000
0.2000 1.1832
0.4000 1.3416
0.6000 1.4832
0.8000 1.6125
1.0000 1.7321
10
例
? ? ?
d 2x dt 2
?
(x2
?
1)
dx dt
?
x
?
0
?? x(0) ? 2; x '(0) ? 0
解: 令 y1=x，y2=y1',
则微分方程变为一阶微分方程组：
? ? ?
y1 y2
' '
? ?
y2 (1 ?
y12 ) y2
?
y1
?? y1 (0) ? 2, y2 (0) ? 0
2、取t0=0, tf=20, 输入命令：
自变 量值
函数 值
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m文件名
ts=[0t，tf]，
t0、tf为自变量 的初值和终
值
函数初 值条件
[t,x]=ode23(‘f', ts, x0级) 2阶3龙格 -库塔公式 [t,x]=ode45(@f, ts, x0) 级4阶5 龙格 -库塔公式
? ? ?
yi? 1 ? yi ? hf y0 ? y( x0 )
( xi
,
yi
)
i ? 0,1,2, ?
,n -1
此即欧拉法(向前欧拉法).
对应有隐式欧拉法
北京邮电大学数学系
? ? ?
yi? 1 ? yi ? hf y0 ? y( x0 )
( xi? 1 ,
yi ? 1 )
3
2、使用数值积分
梯形方法/*trapezoid formula*/
对方程y'=f(x,y), 两边由xi到xi+1积分,并利用梯形公式, 有
? y( xi?1 ) ? y( xi ) ?

xi?1 f ( x, y( x)) dx
xi
?
( xi ? 1
?
xi ) ?
f ( xi ,
y( xi ))
? f ( xi ?1 , 2
y( xi ?1 ))
故有公式
?? ?
yi
因此，研究常微分方程的数值解法十分必要。
对常微分方程
：
? ? ?
y' ? y(x
0
f(x,y) ) ? y0
，其数值解是指由初始点
x 0开始
的若干离散的 xi 处，即对 x0 ? x1 ? x2 ? L ? xn，求出准确值 y(x 1 ),
y( x2 ), L , y( xn ) 的相应近似值 y1 , y2 ,L , yn。
?
1
?
yi
?
h[ 2
f ( xi ,
yi ) ?
f ( xi?1 , yi?1 )]
?? y0 ? y( x0 )
实际应用时，与欧拉公式结合使用
??
y (0 ) i?1
?
yi
?
hf ( xi , yi )
? ??
y(k? i?1
1)
?
yi ?
h 2
[
f
( xi
,
yi
)
?
f
(
xi ? 1
,
y (k ) i?1