微分方程的Matlab求解ppt课件

取t0=0，tf=10，建立程序如下： t0=0;tf=10; [t,y]=ode45('eq3',[t0 tf],[0 0]); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,'-') hold on plot(y(:,1),y(:,2),'*')
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
例4
d2x dt 2
1000(1
x2
)
dx dt
x
0
x(0) 2; x'(0) 0
解: 令 y1=x，y2=y1’
由待解 方程写 成的m-
ts=[t0，tf]， t0、tf为自
函数的 初值
变量的初
文件名 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
z = (-C1+exp(4*t)*C1-2*exp(4*t)+C2+exp(4*t)*C3)*exp(-2*t)
用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s
例 3 求微分方程组的通解.
dx
dt dy
dt dz
dt
2x 4x 4x
3y 5y 4y
3z 3z 2z
解 输命令 ：
s=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't')；
s.x
% 查看结果
s.y
6.大鱼只吃小鱼、小鱼只吃虾米，试建模研究这一捕食系统。在求解你的模型 时也许你会遇到困难，建议对模型中的参数取定几组值，用数值解方法处理， 并研究结果关于参数取值的敏感性。
y1' y2 y3
y2 ' y1 y3 y3' 0.51y1 y2
y1(0) 0, y2 (0) 1, y3(0) 1
解 1、建立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);
x(0) 0, y(0) 0
返回
2. 模型求解
(1) w=20时,建立m-文件如下: function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2); dy(2)=20*(20+15*sin(t)-y(2))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2);
3.从致冰厂购买了一块立方体的冰块，在运输途中发现，第一小时大约融 化了1/4, （1） 求冰块全部融化要多长时间（设气温不变） （2） 如运输时间需要2.5小时，问：运输途中冰块大约会融化掉多少？
4．设药物吸收系数 （k为药物的分解系数），对口服或肌注治疗求体内
药物浓度的峰值（峰浓度）级达峰时间。
s.z x=simple(s.x) %简化结果
y=simple(s.y)
z=simple(s.z)
结 果 为：
x =-(-C1-C2*exp(-3*t)+C2-C3+C3*exp(-3*t))*exp(2*t)
y =(-C1*exp(-4*t)+C1+C2*exp(-4*t)+C2*exp(-3*t)-C2+C3-C3*exp(-3*t))*exp(2*t)
不断修改tf的值,分别取tf=5, 2.5, 3.5,…,至3.15时,狗刚好追上慢跑者.
(2) w=5时
建立m-文件如下: function dy=eq4(t,y) dy=zeros(2,1); dy(1)=5*(10+20*cos(t)-y(1))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2); dy(2)=5*(20+15*sin(t)-y(2))/sqrt ((10+20*cos(t)- y(1))^2+(20+15*sin(t)-y(2))^2);
例 1 求 du 1 u 2 的通解.
dt
解 输入命令：dsolve('Du=1+u^2','t')
结 果：u = tan(t-c)
例 2 求微分方程的特解.
d 2 y 4 dy 29 y 0 dx 2 dx y(0) 0, y' (0) 15
解: 输入命令:
y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 结 果 为 : y =3e-2xsin(5x)
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
2
4
6
8
10
12
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的图形为“+”线.
不断修改tf的值,分别取tf=20, 40, 80,…,可以看出,狗永远追不上慢跑者.
实验题
1. 一个小孩借助长度为a的硬棒，拉或推某玩具.此小孩沿某曲线行走，计 算并画出玩具的轨迹.
研究你画的轨迹的几何形状是由哪些因素确定的？如何确定的？
2. 讨论资金积累、国民收入与人口增长的关系. （1）若国民平均收入x与按人口平均资金积累y成正比，说明仅当总资金积 累的相对增长率k大于人口的相对增长率r时，国民平均收入才是增长的. （2）作出k(x)和r(x)的示意图，分析人口激增会引起什么后果.
0
500
1000
1500
2000
2500
3000
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]);
plot(T,Y(:,1),'-')
3、结果如图
例 5 解微分方程组.
取t0=0，tf=10，建立程序如下： t0=0;tf=10; [t,y]=ode45('eq4',[t0 tf],[0 0]); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,'-') hold on plot(y(:,1),y(:,2),'*')
则微分方程变为一阶微分方程组：
y2
'
y1' y2 1000(1 y12
)
y2
y1
1、建 立my-1(文0)件如2下, y：2 (0) 0
function dy=vdp1000(t,y)
dy=zeros(2,1);
dy(1)=y(2);
2 1.5
1 0.5
0 -0.5
-1 -1.5
-2 -2.5
则X=10+20cost, Y=20+15sint, 狗从(0,0)出发,与导弹追踪问题类似，建立狗 的运动轨迹的参数方程:
dx
dt
w (10 20 cost x)
(10 20 cost x)2 (20 15sin t y)2
dy
dt
w
(20 15sin t y)
(10 20 cost x)2 (20 15sin t y)2
5．正常人身上也有癌细胞，一个癌细胞直径约为10μm，重约0.001μg， （1）当患者被查出患有癌症时，通常直径已有1cm以上（即已增大1000倍）， 由此容易算出癌细胞转入活动期已有30σ天，故如何在早期发现癌症是攻克 癌症的关键之一（2）手术治疗常不能割去所有癌细胞，故有时需进行放射 疗法。射线强度太小无法杀死癌细胞，太强病人身体又吃不消且会使病人免 疫功能下降。一次照射不可能杀死全部癌细胞，请设计一个可行的治疗方案 （医生认为当体内癌细胞数小于100000个时即可凭借体内免疫系统杀灭)。
慢跑者与狗
一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭 圆方程为: x=10+20cost, y=20+5sint. 突然有一只狗攻击他. 这只狗 从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向慢跑者. 分别求出w=20,w=5时狗的运动轨迹.
1. 模型建立
设时刻t慢跑者的坐标为(X(t)，Y(t))，狗的坐标为(x(t)，y(t)).