清华数学实验第七章微分方程与计算机模拟.pptVIP

*/20 微分方程与计算机模拟 常微分方程数值求解方法 追击问题的计算机仿真 有阻力抛射曲线实验 思考题与练习题 ? ? ? ? MATLAB求常微分方程初值问题 数值方法是先创建函数文件,用以描述微分方程右端二元函数,然后用ode23()求出数值解 引例 炮弹在飞行过程中,空气阻力与飞行速度v的平方成正比,如果初始速度v0 ,由牛顿第二定律,得 一阶微分方程主要信息是右端项和初始值: 例7.1 马尔萨斯模型,以1994 年我国人口为12亿为初值，求解常微分方程 N(t)表示人口数量,取人口变化率r =0.015,微分方程 function z=fun1(t,N) z=0.015*N; ode23(fun1,[1994,2020],12) [T,N]=ode23(fun1,[1994,2020],12) 命令窗口 ? 编辑窗口 ? 常微分方程组初值问题 一阶常微分方程组初值问题数值求解方法 [T,y] = ode23( F ,Tspan,y0) 其中, F是函数文件, 表示 微分方程右端函数 Tspan = [t0 Tfinal] —— 求解区域; y0 —— 初始条件 注: 函数F(t,y) 必须返回列向量. 数值解 y 的每一行对应于列向量T中的每一行数据 捕食者与被捕食者问题 海岛上有狐狸和野兔,当野兔数量增多时，狐狸捕食野兔导致狐群数量增长;大量兔子被捕食使狐群进入饥饿状态其数量下降;狐群数量下降导致兔子被捕食机会减少,兔群数量回升。微分方程模型如下 计算 x(t)，y(t) 当t∈[0，20]时的数据。绘图并分析捕食者和被捕食者的数量变化规律。 x(0)= 100 y(0)=20 创建MATLAB的函数文件 function z=fox(t,y) z(1,:)=y(1)-0.015*y(1).*y(2); z(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) 求微分方程数值解并绘解函数图形 ------兔子数量; ------狐狸数量 兔-狐数量 变化相位图 例7.3 “蝴蝶效应”来源于洛伦兹一次讲演。模型如下 求微分方程数值解, 绘出解函数在Y-X平面投影曲线 取 ?=8/3，?=10，?=28。 x(0)=0，y(0)=0，z(0)=0.01。 t∈[0，80]， 记向量 [y1，y2，y3] = [x，y，z]，创建函数文件 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); figure(1),plot(y(:,2),y(:,1)) figure(2),comet3(y(:,1),y(:,2),y(:,3)) 例7.7 追击问题模拟。设系统中有动点Q和动点P，点Q从坐标原点出发以速度V=1(米/秒)沿y轴正向匀速直线运动，点P从坐标原点右侧100米处与Q点同时出发，以2V速度紧盯Q点追赶。60秒后P能否追上Q 。 时间步长法——模拟随时间变化的系统状态 计算机模拟图 ? 时间以步长dt 向前推进时,系统中两个动点在各个时刻的速度、位移、位置和两点间距离 在平面坐标系中, 初始时刻点Q的坐标(0, 0), 点P的坐标(100, 0)。在时刻 tk , 点Q 以均匀速度v=1(m/min)沿Y轴正向运动,而点P以2v的速度追赶Q。 x O 100 60 y 时刻 tk ,点Q坐标为: (0, vk),点P的坐标为: (xk, yk) 追击方向 位置 function 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; end x=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位置向量 %计算追赶方向 %追赶方向单位