9传染病模型与微分方程数值解

四阶龙格-库塔公式
6
4,数值公式的精度 , 当一个数值公式的截断误差可表示为o(hk)时 (k为正整数,h为步长),称它是一个k阶公式 阶公式. 阶公式 k越大,则数值公式的精度越高. 欧拉法是一阶公式,改进的欧拉法是二阶公式.
线性多步法有四阶阿达姆斯外插公式和内插公式.
北京邮电大学数学系
7
(三)用Matlab软件求常微分方程的数值解 三用 软件求常微分方程的数值解 [t,x]= solver( 'f', ts, x0, options) , )
h yi+1 = yi + [ f ( xi , yi ) + f ( xi+1 , yi+1 )] 2 y0 = y( x0 )
故有公式
实际应用时,与欧拉公式结合使用
0 y i(+1) = y i + hf ( xi , y i ) h ( k +1) k y i +1 = y i + [ f ( xi , y i ) + f ( xi +1 , y i(+1) )] k = 0,1,2, 2
第五章 微分方程模型
常微分方程的数值解
5.1 传染病模型
北京邮电大学数学系
1
常微分方程的数值解及实验
(一)常微分方程数值解 一 常微分方程数值解 高数中微分方程解法在实际中基本不会直接使用 在生产和科研中所处理的微分方程往往很复杂且大多 得不出一般解.而在实际上对初值问题,一般是要求 得到解在若干个点上满足规定精确度的近似值 解在若干个点上满足规定精确度的近似值,或者 解在若干个点上满足规定精确度的近似值 得到一个满足精确度要求的便于计算的表达式 满足精确度要求的便于计算的表达式. 满足精确度要求的便于计算的表达式 因此,研究常微分方程的数值解法十分必要. 因此,研究常微分方程的数值解法十分必要
y i +1 = y i + hf ( xi , y i ) i = 0,1,2, , n - 1 y 0 = y ( x0 )
此即欧拉法 向前欧拉法). 欧拉法(向前欧拉法 欧拉法 向前欧拉法
对应有隐式欧拉法 对应有隐式欧拉法
yi +1 = yi + hf ( xi +1 , yi +1 ) y0 = y ( x0 )
10
例
d 2x dx 2 2 + ( x 1) + x = 0 dt dt x(0) = 2; x '(0) = 0
解: 令 y1=x,y2=y1', 则微分方程变为一阶微分方程组:
[t,y]=ode45(@dfun2,[0,20],[2,0]); [t,y] 1,建立 文件 文件dfun2.m如下: plot(t,y(:,1),'r-',t,y(:,2),'b.-'); 如下: ,建立m-文件 如下 function dx=dfun2(t,y) hold on dx=[y(2); plot(y(:,1),y(:,2),'co'); (1-y(1)^2)*y(2)-y(1)]; hold off legend('t ~ x', 't ~ x`', 'x ~ y');
(二)建立数值解法的一些途径 二 建立数值解法的一些途径
设 x i +1 xi = h, i = 0,1,2, n 1, 可用以下离散化方法求解微分方程: y' = f(x, y) y(x 0 ) = y 0
1,用差商代替导数 , 若步长h较小,则有 故有公式:
y ( x + h) y ( x ) y ' ( x) ≈ h
3,结果如图
d 2x dx 2 2 + ( x 1) + x = 0 dt dt x(0) = 2; x '(0) = 0
y1 ' = y2 y2 ' = (1 y12 ) y2 y1 y (0) = 2, y (0) = 0 2 1
北京邮电大学数学系
12
动态模型
描述对象特征随时间 空间 的演变过程 描述对象特征随时间(空间 空间)的演变过程 分析对象特征的变化规律 预报对象特征的未来性态 研究控制对象特征的手段 微分 方程 根据建模目的和问题分析作出简化假设 建模 按照内在规律或用类比法建立微分方程
北京邮电大学数学系 15
模型一 模型假设: 模型假设:
(1)一人得病后,久治不愈,人在传染期内不会死亡. 病人) 设已感染人数 (病人 x(t), 假设是连续可微函数 病人 (2) 假设每个病人每天有效接触(足以使人致病)人数为λ
建模
x(t + t ) x(t ) = λ x(t )t
dx = λx dt x (0 ) = x 0
3级2阶龙格 库塔公式 级 阶龙格 阶龙格-库塔公式 5级4阶龙格 库塔公式 级 阶龙格 阶龙格-库塔公式
用于设定误差限 误差限(缺省时设定相对误差10-3, 绝对误差10-6), 误差限 命令为:options=odeset('reltol',rt,'abstol',at), rt, at:分别为设定的相对误差 绝对误 相对误差和绝对误 相对误差 绝对误差.
x (t ) = x0 e λ t
t →∞ x →+∞ ?
16
北京邮电大学数学系
举个实例 最初只有1个病人,1个病人一天可传染1个人
x=0:0.1:10; y=exp(x); plot(x,y,'b-');
exp(10) = 22026
北京邮电大学数学系
17
t →∞ x →+∞
模型缺点
问题: 问题:随着时间的推移,病人的数目将无限增加, 这一点与实际情况不符. 原因: 原因: 当不考虑传染病期间的出生,死亡和迁移时, 一个地区的总人数可视为常数, 因此 λ应为时间t的函数. 在传染病流行初期, λ较大, 随着病人的增多,健康人数减少, 被传染的机会也减少,于是λ将变小. 模型修改的关键: 模型修改的关键:λ的变化规律 若有效接触的是病人, 若有效接触的是病人,则不能使病人数增加 必须区分已感染者(病人 和未感染者 健康人) 必须区分已感染者 病人)和未感染者 健康人 病人 和未感染者(健康人
北京邮电大学数学系 8
注意: 注意 1,在解n个未知函数的方程组时,x0和x均为n维向 量,m-文件中的待解方程组应以x的分量形式写成.
2,使用Matlab软件求数值解时,高阶微分方程必 须等价地变换成一阶微分方程组.
北京邮电大学数学系
9
微分方程求解实例
设取步长 ,从 到 用四阶龙格-库塔方法
求解初值问题 建立m-文件dfun1.m如下 function dx=dfun1(x,y) dx=y-2*x/y; 输入命令
自变 量值 函数 值 ode45 ode23 ode113 ode15s ode23s 由待解 方程写 成的m成的 文件名 ts=[t0,tf], , t0,tf为自变 , 量的初值和 终值 函数初 值条件
[t,x]=ode23('f', ts, x0) [t,x]=ode45(@f, ts, x0)
5
3,使用泰勒公式 , 以此方法为基础,有龙格 库塔法 线性多步法 龙格-库塔法 线性多步法等方法 龙格 库塔法,线性多步法 库塔三阶方法 库塔
h yn+1 = yn + 6 (K1 + 4K2 + K3 ), K1 = f (xn , yn ), K = f (x + h , y + h K ), n n 1 2 2 2 K = f (x + h, y hK + 2hK ). n n 1 2 3
y′( x1 ) ≈ y ( x 2 ) y ( x0 ) 2h
y x1 ) (
y ( x 2 ) ≈ y ( x 0 ) + 2 h f ( x1 , y ( x1 ))
x0
y i +1 = yi 1 + 2h f ( x i , y i ) i = 1, ... , ni (t + t ) i (t )] = [ Ni (t ) λ s (t )]t di di = λ si = λ i (1 i ) dt dt s (t ) + i (t ) = 1 i ( 0 ) = i0
北京邮电大学数学系
19
di = λ i (1 i ) dt i ( 0 ) = i0
感染疾病的人数与哪些因素有关? 感染疾病的人数与哪些因素有关 描述传染病的传播过程 分析受感染人数的变化规律 预报传染病高潮到来的时刻 预防传染病蔓延的手段
北京邮电大学数学系 14
问题分析 不同类型传染病的传播过程有不同的特点, 故不从 医学的角度对各种传染病的传播过程一一进行分析, 而是按一般的传播机理建立模型. 由于传染病在传播的过程涉及因素较多,在分 析问题的过程中,不可能通过一次假设建立完善的 数学模型. . 思路: 思路 先做出最简单的假设,对得出的结果进行分析, 针对结果中的不合理之处 结果中的不合理之处,逐步修改假设, 结果中的不合理之处 最终得出较好的模型.
北京邮电大学数学系 11
y1 ' = y2 y2 ' = (1 y12 ) y2 y1 y (0) = 2, y (0) = 0 2 1
2,取t0=0, tf=20, 输入命令:
[t,x]=ode45(@dfun2,[0,20],[2,0]) plot(t,x(:,1),'r-',t,x(:,2),'b.-'); hold on plot(x(:,1),x(:,2),'co'); hold off legend('t ~ x','t ~ x`','x ~ y');