第1讲 微分方程模型

1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6
解
1、建立m-文件rigid.m如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
-0.8 -1 0 2 4 6 8 10 12
例3.求解下列初值问题：
解： （1）创建函数文件：xprim2.m
x1 ' x1 0.1x1 x2 0.01t , x1 (0) 30, x2 ' x2 0.02 x1 x2 0.04t , x2 (0) 20.
function y=xprim2(t,x) y=[x(1)-0.1*x(1)*x(2)+0.01*t; -x(2)+0.02*x(1)*x(2)+0.04*t]; （2）创建求解函数文件：solplot.m（脚本文件） ，文件中 调用matlab的函数ode45求方程的数值解，并画出解曲线及相图。
b. 微分方程的求解方法---解析解与数值解
（2）数值解法
数学软件：matlab, mathematica, maple Matlab：ode45（非刚性方程） ode23（刚性方程） Mathematica：ode45（非刚性方程） ode23（刚性方程）
用Matlab软件求常微分方程的数值解
b. 微分方程的求解方法---解析解与数值解
（1）解析解法---初等积分法
一阶微分方程：可分离变量, 齐次方程, 线性方程等. '' ' 二阶齐次微分方程: ax bx cx 0
b b2 4ac b b 2 4ac , r2 解特征方程 r 1 2a 2a
b 2 4ac 0, b 4ac 0,
这个结论对 于(4-1)也是 成立的.
常微分方程组的平衡点及其稳定性
定义：设
dx dt f ( x, y ), dy g ( x, y ). dt
f ( x, y ) 0, g ( x, y ) 0.
(3)
代数方程组
的实根x = x0, y = y0称为方程(3)的平衡点, 记作P0 (x0, y0). 它也是方程(3)的解.
如果
t
lim x(t ) x0 ,
t
lim y(t ) y0 ,
则称平衡点P0是稳定的.
下面给出判别平衡点P0是否稳定的判别准则. 定理：设 f f x ( P0 ) y ( P0 ) f q det g , p ( P0 ) ( P0 ) , f g y x x ( P0 ) y ( P0 ) 则 （1）当p＞0且q＞0时, 平衡点P0是稳定的； （2）当p＜0或q＜0时, 平衡点P0是不稳定的.
用Matlab求微分方程的解析解
例1：求微分方程 解：输入命令: dsolve(‘Du=1+u^2’,’t’) 输出结果: u=tan(t+c1)
du 1 u2 dt
的通解。
例2
求微分方程的特解.
d 2 y dy 2 4 29y 0 dx dx y (0) 0, y ' (0) 15
-0.5 -1 -1.5 -2 -2.5 0
y1 ' y2 y2 ' 1000 1 y12 ) y2 y1 ( y (0) 2, y (0) 0 2 1
2 1.5 1 0.5 0
500
1000
1500
2000
2500
3000
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-') 3、结果如图
2
x c1e r1t c2e r2t b x c1e c2te , r1 r2 2a x e t (c1 cos t c2 sin t )
r1t r1t
b 2 4ac 0,
b , 2a
4ac b 2 2a
用Matlab求微分方程的解析解
[t，x]=solver（’f’,ts,x0,options）
自变 量值 函数 值
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m文件名
ts=[t0，tf]， 函数的 初值 t0、tf为自 变量的初 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法 用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
动态模型
• 描述对象特征随时间的演变过程 • 分析对象特征的变化规律 • 预报对象特征的未来性态 • 研究控制对象特征的手段
常用动态模型
• 微分方程模型（连续） • 差分方程模型（离散） • 时间序列模型（随机）
第一讲
微分方程模型பைடு நூலகம்
1.0 常微分方程基本理论与建模方法
1.1 传染病模型
1.2 经济增长模型 1.3 药物在体内的分布与排除 1.4 人口预测和控制 1.5 烟雾的扩散与消失
例2
解: 令 y1=x，y2=y1’
d 2x dx 2 1000 1 x 2 ) x 0 ( dt dt x(0) 2; x' (0) 0
则微分方程变为一阶微分方程组：
1、建立m-文件vdp1000.m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
来代替.
dx f ( x0 )( x x0 ) dt
(2)
易知 x0也是方程(2)的平衡点. (2)的通解为
x(t ) Ce
f ( x0 ) t
关于x0是否稳定有以下结论：
x0 ,
① 若 f ( x0 ) 0, 则x0是稳定的；
② 若 f ( x0 ) 0, 则x0是不稳定的.
动态模型
Dynamical Models
南京农业大学 数学系 2011年8月
数学建模的一般步骤
1. 了解问题的实际背景，明确建模目的，收集掌握 必要的数据资料。 2. 通过对资料的分析计算，找出起主要作用的因素， 经必要的精练、简化，提出若干符合客观实际的 假设。 3. 在所作假设的基础上，利用适当的数学工具去刻 画各变量之间的关系，即建立模型。 4. 模型求解（包括解方程、图解、逻辑推理、定理 证明等）。 5. 模型的分析与检验。
dx f (t , x), dt x(t0 ) x0 , • 一般高阶线性常微分方程具有形式 dy0 n n 1 y1 , d x dx d x f ( t ; x , ,..., n1 ) dt n dt dt dt dy1 y , i 2 d x dt 令 yi ( i 0,1,...,n) i dt .... dyn1 f ( t ; y , y ,..., y ) • 通解 0 1 n 1 dt • 初值问题（IVP）与特解
例1. 求解下列初值问题： x ' x 2 , x (0) 1. 解： （1）创建函数文件：xprim1.m
function y=xprim1(t,x) y=-x.^2; （2）创建求解函数文件：solplot.m（脚本文件） ，文件中 调用matlab的函数ode45求方程的数值解，并画出解曲线。 [t, x]=ode45(‘xprim1’, [0,1],1); plot(t,x,’-’,t,x,’o’); xlabel(‘Time: t0=0, tend=1’); ylabel(‘x values: x(0)=1’);
定性分析方法（几何方法）：不求方程的解，而直 接分析解的（几何）性质。 定义：设
dx f ( x) dt (1)
称代数方程 f (x)=0 的实根x = x0为方程的平衡点 (或奇点). 它也是方程的解. 如果
lim x(t ) x0
t
则称平衡点x0是稳定的.
稳定性判别方法 由于 f ( x) f ( x )( x x ), 在讨论方程(1)的 0 0 稳定性时，可用其线性近似
[t, x]=ode45(‘xprim2’, [0,20],[30;20])
plot(t,x); xlabel(‘Time: t0=0, tend=20’); ylabel(‘x values: x1(0)=30,x2(0)=20’);
例 4 解微分方程组. y1 ' y2 y3 y2 ' y1 y3 y3 ' 0.51y1 y2 y1 (0) 0, y2 (0) 1, y3 (0) 1
解:
输入命令 ：
[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't')； x=simple(x) % 将x化简 y=simple(y) z=simple(z)
输出结 果 ：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
命令：
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
记号: 在表达微分方程时，用字母D表示求微分，D2、 D3等表示求高阶微分. 任何D后所跟的字母为因变量， 自变量可以指定或由系统规则选定为确省.
d2 y 例如，微分方程有 2 0 dx
应表达为：D2y=0.
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图 图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
C. 常微分方程的定性方法----相图与稳定性分析
解: 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 输出结果 : y =3e-2xsin（5x）
例3
求微分方程组的通解. dx dt 2 x 3 y 3z dy 4 x 5 y 3z dt dz 4 x 4 y 2 z dt
•定义：
常微分方程解的存在唯一性
定理：对初值问题（IVP）
dx f (t , x) dt x(t0 ) x0
(1) (2)
设函数f (t , x)在矩形域R {(t , x) :| t t0 | a,| x x0 | b} 上连续,且关于x满足Lipschitz条件. 即存在常数L 0, 使得 | f (t , x1 ) f (t , x2 ) | L | x1 x2 |, (t , x1 ), (t , x2 ) R. 则方程(1)在区间 | t t0 | h上存在满足(2)的解x (t ) b 且解唯一.其中h min(a, ), M max | f ( x) | . R m
1.6 万有引力定律的发现
1.0 常微分方程基本理论 与建模方法
1.0.1 常微分方程的基本问题
a. 常微分方程解的存在唯一性
b. 常微分方程的求解方法----初等积分法与数值方法
c. 常微分方程的定性方法----相图与稳定性分析
d. 常微分方程建模----两个简单实例
a. 常微分方程基本概念及解的存在唯一性