《清华数学实验第七章微分方程与计算机模拟》由会员分享，可在线阅读，更多相关《清华数学实验第七章微分方程与计算机模拟（20页珍藏版）》请在金锄头文库上搜索。
1、微分方程与计算机模拟,常微分方程数值求解方法 追击问题的计算机仿真 有阻力抛射曲线实验 思考题与练习题, ,MATLAB求常微分方程初值问题,数值方法是先创建函数文件,用以描述微分方程右端二元函数,然后用ode23()求出数值解,引例 炮弹在飞行过程中,空气阻力与飞行速度v的平方成正比,如果初始速度v0 ,由牛顿第二定律,得,一阶微分方程主要信息是右端项和初始值:,例7.1 马尔萨斯模型,以1994 年我国人口为12亿为初值，求解常微分方程,N(t)表示人口数量,取人口变化率r =0.015,微分方程,function z=fun1(t,N) z=0.015*N;,ode23(fun1,199
2、4,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中的每一行数据,捕食者与被捕食者问题,海岛上有狐狸和野兔,当野兔数量增多时，狐狸捕食野兔导致狐群数量增长;大量兔子被捕食使狐群进入饥饿状态其数量下降;狐群数量下降导致兔子被捕食机会减少,兔群数量回
3、升。微分方程模型如下,计算 x(t)，y(t) 当t0，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),求微分方程数值解并绘解函数图形,-兔子数量; -狐狸数量,
4、兔-狐数量 变化相位图,例7.3 “蝴蝶效应”来源于洛伦兹一次讲演。模型如下,求微分方程数值解, 绘出解函数在Y-X平面投影曲线,取 =8/3，=10，=28。,x(0)=0，y(0)=0，z(0)=0.01。,t0，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,
5、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 以
6、均匀速度v=1(m/min)沿Y轴正向运动,而点P以2v的速度追赶Q。,时刻 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 (米),%设置初值 %设置
7、追赶方向 %计算P位置向量 %计算Q位置向量 %计算追赶方向 %追赶方向单位化 %提取追击曲线坐标 %绘追击曲线,追击问题动态模拟程序,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) end,抛射曲线实验，假设阻力与速度成正比。在微分方程中增加阻力项,符号计
8、算方法,syms t v g alfa k x=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),2008年贺岁片集结号展现出视听震撼的战争场面，电影中92式山炮，炮弹初速: 198米/秒,最大射程:2788米 利用实验程序确定阻力系数 k,function Xmax=mlab72(k) alfa=pi/4; v=198;g=9.8; t=0;dt=.1; x=0;y=0; while
9、 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,xk;y=y,yk; end Xmax=xk; plot(x,y,ro),实验数据:,k=0.02,k=0.015,例7.10 设圆柱直径等于圆球半径，圆柱面紧贴球的中轴线穿过球面时留下切痕。 模拟圆柱面切割球面过程,并绘出切割剩余的空间立体,圆柱面方程,球面方程,圆柱面用极坐标表示为,function viviana() theta=(-180:5:180)*pi/180; fai=(
10、-90:5:90)*pi/180; X=cos(theta)*cos(fai); Y=cos(theta)*sin(fai); Z=sin(theta)*ones(size(fai); colormap(0 0 1) mesh(X,Y,Z),axis off hold on,theta=(360:-10:0)*pi/180; x=.5*cos(theta)+.5; y=.5*sin(theta); U=x;x;V=y;y; E=ones(size(x); W=-E;1.2*E; for p=10:fix(360/dt)+1 II=1:p; u=U(:,II); v=V(:,II); w=W(:,II); mesh(u,v,w),pause(.5) end,思考题与练习题,蛇形曲线的微分方程 右端函数 在平面区域内任意一点(x,y)处的值确定解曲线的切 线斜率.利用quiver(x,y,px,py)绘平面向量场,3. 对于有阻力的抛射体的抛射曲线参数方程,推导飞行时间 T, 以及飞行距离 X 的数学表达式,2. 在追击问题的模拟程序中,为了计算表示追击方向的单位向量e, 使用了命令norm()。利用help norm 了解norm的主要功能以及norm()更多的使用方法,
