微分方程与计算机模拟课件.ppt

1、微分方程与计算机模拟微分方程与计算机模拟 常微分方程数值求解方法常微分方程数值求解方法追击问题的计算机仿真追击问题的计算机仿真有阻力抛射曲线实验有阻力抛射曲线实验思考题与练习题思考题与练习题2/19例例7.1 .1 以以 1994 年我国人口为年我国人口为12亿为初值，求解马尔亿为初值，求解马尔萨斯模型常微分方程萨斯模型常微分方程 N(t)表示人口数量表示人口数量, ,取人口变化率取人口变化率r = =0.015, ,微分方程微分方程NdtdN015. 0 function z=fun1(t,N)z=0.015*N; ode23(fun1,1994,2020,12)T,N=ode23(fun1

2、,1994,2020,12)199019952000200520102015202012141618命令窗口命令窗口 编辑窗口编辑窗口 12)1994( N3/19例例2. Logistic模型模型 )1(uudtdu 0)0(uu 创建微分方程右端函数创建微分方程右端函数:function z=fun2(t,u)z=u.*(1-u);ode23(fun2,0,6,1.8)ode23(fun2,0,6,0.2)在命令窗口求数值解在命令窗口求数值解4/19x,y=meshgrid(0:.25:6,0:.05:2);k=y.*(1-y);d=sqrt(1+k.2);px=1./d;py=k./d;

3、quiver(x,y,px,py),hold onu=dsolve(Du=u*(1-u),u(0)=.2); v=dsolve(Dv=v*(1-v),v(0)=1.8);ezplot(u,0,6)ezplot(v,0,6)根据微分方程右端函数根据微分方程右端函数 f(x，y)= y(1 y)，区域区域D = (x, y) | 0 x 6, 0 y 2内未知函数的导数值，确定解函数曲线的切线对应单内未知函数的导数值，确定解函数曲线的切线对应单位向量，绘制向量场。位向量，绘制向量场。5/19常微分方程组初值问题常微分方程组初值问题一阶常微分方程组初值问题数值求解方法一阶常微分方程组初值问题数值求解

4、方法T,y = ode23( F ,Tspan,y0) 其中其中, F是函数文件是函数文件, 表示表示 微分方程右端函数微分方程右端函数Tspan = t0 Tfinal 求解区域求解区域; y0 初始条件初始条件注注: 函数函数F(t,y) 必须返回列向量必须返回列向量. 数值解数值解 y 的每一行对应于列向量的每一行对应于列向量T中的每一行数据中的每一行数据 000)(),(ytyttytfdtyd6/19捕食者与被捕食者问题捕食者与被捕食者问题 海岛上有狐狸和野兔海岛上有狐狸和野兔, ,当野兔数量增多时，狐狸捕食当野兔数量增多时，狐狸捕食野兔导致狐群数量增长野兔导致狐群数量增长; ;大量

5、兔子被捕食使狐群进入大量兔子被捕食使狐群进入饥饿状态其数量下降饥饿状态其数量下降; ;狐群数量下降导致兔子被捕食狐群数量下降导致兔子被捕食机会减少机会减少, ,兔群数量回升。微分方程模型如下兔群数量回升。微分方程模型如下 xyydtdyxyxdtdx01. 0015. 0计算计算 x(t)，y(t) 当当t0，20时的数据。绘图并分时的数据。绘图并分析捕食者和被捕食者的数量变化规律。析捕食者和被捕食者的数量变化规律。x(0)= 100y(0)=20 7/19创建创建MATLAB的函数文件的函数文件function z=fox(t,y)z(1,:)=y(1)-0.015*y(1).*y(2);z

6、(2,:)=-y(2)+0.01*y(1).*y(2);Y0=100,20;t,Y=ode23(fox,0,20,Y0);x=Y(:,1);y=Y(:,2);figure(1),plot(t,x,b,t,y,r)figure(2),plot(x,y) 求微分方程数值解并绘解函数图形求微分方程数值解并绘解函数图形8/19-兔子数量; -狐狸数量兔-狐数量变化相位图9/19“蝴蝶效应蝴蝶效应”来源于洛伦兹一次讲演来源于洛伦兹一次讲演。模型如下模型如下 zyxydtdzzydtdyyzxdtdx )(求微分方程数值解求微分方程数值解, , 并绘出解函数在并绘出解函数在Y-X平面的投影平面的投影曲线曲

7、线 取取 =8/3， =10， =28。x(0)=0，y(0)=0，z(0)=0.01。t0，80，-20-1001020020406010/19记向量记向量 y1，y2，y3 = x，y，z，创建创建MATLAB函数函数文件如下文件如下 function z=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3);z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用用MATLAB命令求解并绘出命令求解并绘出Y-X平面的投影图平面的投影图 y0=0;0;0.01;x,y=ode23(flo,0, 80,y0);figur

8、e(1),plot(y(:,2),y(:,1) figure(2),comet3(y(:,1),y(:,2),y(:,3)11/19例例7.7 追击问题模拟追击问题模拟。设系统中有动点设系统中有动点Q和动点和动点P，点，点Q从坐标原点出发以速度从坐标原点出发以速度V=1(米米/ /秒秒)沿沿y轴正向匀速直轴正向匀速直线运动，点线运动，点P从坐标原点右侧从坐标原点右侧100米处与米处与Q点同时出发，点同时出发，以以2V速度紧盯速度紧盯Q点追赶。点追赶。60秒后秒后P能否追上能否追上Q 。时间步长法时间步长法模拟随时间变化的系统状态模拟随时间变化的系统状态计算机模拟图计算机模拟图 时间以步长时间以

9、步长dt 向前推向前推进时进时, ,系统中两个动系统中两个动点在各个时刻的速度、点在各个时刻的速度、位移、位移、位置位置和两点间和两点间距离距离12/19在平面坐标系中在平面坐标系中, 初始时刻点初始时刻点Q的坐标的坐标(0, 0), 点点P的坐标的坐标(100, 0)。在时刻在时刻 tk , 点点Q 以均匀速度以均匀速度v=1(m/min)沿沿Y轴轴正向运动正向运动,而点而点P以以2v的速度追赶的速度追赶Q。xO10060 y时刻时刻 tk , ,点点Q坐标为坐标为: : (0, vk), ,点点P的坐标为的坐标为: (xk, yk),0)(122kkkkkkkyvxyvxe kkkkkey

10、xyx2,11 追击方向追击方向位置位置13/19function d = chase()P=100,0;Pk=P;Q=0,0;e=-1,0;for k=1:60 Pk=Pk+2*e;P=P;Pk; Qk=0,k;Q=Q;Qk; e=Qk-Pk; d=norm(e);e=e/d;endx=P(:,1);y=P(:,2);u=Q(:,1);v=Q(:,2);plot(u,v,o,x,y,r*) 60秒后秒后P点与点与Q点的距离点的距离 d = 7.0619 (米米)%设置初值设置初值%设置追赶方向设置追赶方向%计算计算P位置向量位置向量%计算计算Q位置向量位置向量%计算追赶方向计算追赶方向%追赶

11、方向单位化追赶方向单位化%提取追击曲线坐标提取追击曲线坐标%绘追击曲线绘追击曲线14/19追击问题动态模拟程序追击问题动态模拟程序function d = chase()Pk=100,0;P=Pk;Q=0,0;e=-1,0;for k=1:60 Pk=Pk+2*e;P=P;Pk; Qk=0,k;Q=Q;Qk; e=Qk-Pk; d=norm(e);e=e/d; x=P(:,1);y=P(:,2); u=Q(:,1);v=Q(:,2); plot(u,v,o,x,y,r*,0,60,og),pause(.5)end020406080100010203040506015/19抛射曲线实验抛射曲线实

12、验，假设阻力与速度成正比。在微分方，假设阻力与速度成正比。在微分方程中增加阻力项程中增加阻力项 )()()()(tykgtytxktx符号计算方法符号计算方法 sin)0(, 0)0(cos)0(, 0)0(00vyyvxx syms t v g alfa kx=dsolve(D2x=-k*Dx,x(0)=0,Dx(0)=v*cos(alfa);y=dsolve(D2y=-g-k*Dy,y(0)=0,Dy(0)=v*sin(alfa);X=taylor(x,3,t),Y=simplify(taylor(y,3,t) sin2121sin)(cos21cos)(222vtkgttvtyvkttv

13、tx16/192008年贺岁片年贺岁片集结号集结号展现出视听震撼的战争场展现出视听震撼的战争场面，电影中面，电影中92式山炮式山炮，炮弹初速炮弹初速:198米米/秒秒,最大射程最大射程:2788米米利用实验程序确定阻力系数利用实验程序确定阻力系数 kfunction Xmax=mlab72(k)alfa=pi/4; v=198;g=9.8;t=0;dt=.1;x=0;y=0;while y=0; t=t+dt; xk=v*cos(alfa)*t-1/2*v*cos(alfa)*k*t2; yk=v*sin(alfa)*t-1/2*g*t2-1/2*t2*v*sin(alfa)*k; x=x,x

14、k;y=y,yk;endXmax=xk;plot(x,y,ro)17/1905001000150020002500-2000200400600800050010001500200025003000-20002004006008001000实验数据实验数据:k 0.1 0.01 0.02 0.015Xmax 677.35 3073.15 2433.66 2719.33k=0.02k=0.01518/19思考题与练习题思考题与练习题1.1. 蛇形曲线的微分方程蛇形曲线的微分方程 右端函数右端函数在平面区域内任意一点在平面区域内任意一点(x,y)处的值确定解曲线的切处的值确定解曲线的切线斜率线斜率.

15、 .利用利用quiver(x,y,px,py)绘平面向量场绘平面向量场yxy2112 -6-4-20246-1.5-1-0.500.5119/193. 对于有阻力的抛射体的抛射曲线参数方程对于有阻力的抛射体的抛射曲线参数方程 sin2121sin)(cos21cos)(222vtkgttvtyvkttvtx推导飞行时间推导飞行时间 T, 以及飞行距离以及飞行距离 X 的数学表达式的数学表达式2. 在追击问题的模拟程序中在追击问题的模拟程序中,为了计算表示追击方为了计算表示追击方向的单位向量向的单位向量e, 使用了命令使用了命令norm()。利用。利用help norm 了解了解norm的主要功

16、能以及的主要功能以及norm()更多的更多的使用方法使用方法20/19附：附：蝴蝶效应蝴蝶效应 蝴蝶效应这一说法来源，是气象学家洛仑兹制作了一个电脑程蝴蝶效应这一说法来源，是气象学家洛仑兹制作了一个电脑程序，模拟气候的变化，并用图像来表示。最后他发现，图像是序，模拟气候的变化，并用图像来表示。最后他发现，图像是很敏感地依赖于初值，而且十分像一只蝴蝶张开的双翅，因而很敏感地依赖于初值，而且十分像一只蝴蝶张开的双翅，因而他形象的将这一图形以他形象的将这一图形以“蝴蝶扇动翅膀蝴蝶扇动翅膀”的方式进行阐释。蝴的方式进行阐释。蝴蝶效应是指对初始条件敏感性的一种依赖现象。输入端微小的蝶效应是指对初始条件敏

17、感性的一种依赖现象。输入端微小的差别会迅速放大到输出端。蝴蝶效应在经济生活中比比皆是：差别会迅速放大到输出端。蝴蝶效应在经济生活中比比皆是：中国宣布发射导弹，港台中国宣布发射导弹，港台100亿美元流向美国。今天的蝴蝶效亿美元流向美国。今天的蝴蝶效应应”或者或者“广义的蝴蝶效应广义的蝴蝶效应”已不限于当初洛仑兹的蝴蝶效应已不限于当初洛仑兹的蝴蝶效应仅对天气预报而言，而是一切复杂系统对初值极为敏感性的代仅对天气预报而言，而是一切复杂系统对初值极为敏感性的代名词或同义语。名词或同义语。 核心理念：看似微不足道的细小变化，却能以某种方式对社会核心理念：看似微不足道的细小变化，却能以某种方式对社会产生微妙的影响，甚至影响整个社会系统的正常运行。关注细产生微妙的影响，甚至影响整个社会系统的正常运行。关注细节，防微杜渐，注重关联，控制全局。节，防微杜渐，注重关联，控制全局。