第七章微分方程与计算机模拟注记

第七章微分方程与计算机模拟PPT注记（2010-3-24）只有三个内容，有一定难度和深度。三个内容是：常微分方程初值问题求数据解和蝴蝶效应（洛仑兹模型）的动态仿真、追击曲线动态仿真、有阻力的抛射曲线实验中阻力系数（电影集结号中的火炮数据）一、一阶常微分方程MATLAB求常微分方程初值问题00)(),(yxyyxfy数值解方法是首先创建一个函数文件，用以描述一阶微分方程右端的二元函数，然后用ode23()求解对应的微分方程初值问题。由于命令行较少，可以在命令窗口实现第一条命令ode23('fun1',[1994,2020],12)执行后，在图形窗口显示解曲线；第二条命令[T,N]=ode23('fun1',[1994,2020],12)执行后，在命令窗口显示数值解的数据与练习题有关的两个例题马尔萨斯模型：NdtdN15.0，]3,0[]3,0[),(yx00.511.522.533.5-0.500.511.522.533.5clear[x,y]=meshgrid(0:0.2:3);k=0.15*y;d=sqrt(1+k.^2);px=1./d;py=k./d;quiver(x,y,px,py)axisimage蛇形曲线模型：yxy2112[x,y]=meshgrid(-5:.5:5,-1:.1:1);k=1./(1+x.^2)-2*y.^2;d=sqrt(1+k.^2);px=1./d;py=k./d;quiver(x,y,px,py)-6-4-20246-1.5-1-0.500.51二、洛仑兹模型的动态演示蝴蝶效应这一说法来源，是气象学家洛仑兹制作了一个电脑程序，模拟气候的变化，并用图像来表示。最后他发现，图像是很敏感地依赖于初值，而且十分像一只蝴蝶张开的双翅，因而他形象的将这一图形以“蝴蝶扇动翅膀”的方式进行阐释。蝴蝶效应是指对初始条件敏感性的一种依赖现象。输入端微小的差别会迅速放大到输出端。蝴蝶效应在经济生活中比比皆是：中国宣布发射导弹，港台100亿美元流向美国。今天的蝴蝶效应”或者“广义的蝴蝶效应”已不限于当初洛仑兹的蝴蝶效应仅对天气预报而言，而是一切复杂系统对初值极为敏感性的代名词或同义语。核心理念：看似微不足道的细小变化，却能以某种方式对社会产生微妙的影响，甚至影响整个社会系统的正常运行。关注细节，防微杜渐，注重关联，控制全局。comet()平面曲线慧星绘图，用于模拟平面曲线形成的动态过程；comet3()是三维曲线慧星图绘图命令，用于模拟三维曲线形成的动态过程。这两条命令的弱点是对绘图数据量需求很大，数据量不够大时，产生的动态过程时间太短，没有动态过程的效果。222222((),()),(0,)((),())()(1)()4(2)(0)100,(0)0()()(),(),4ttttttxxxxxxxttPxtytQtPQxttytxxtytytxyxyytytxxtyyxttxyytytyxxxyttytyxxxxy三、追击曲线程序02040608010002040601．静态追击曲线图形绘制程序functiond=chase()Pk=[100,0];P=Pk;Q=[0,0];e=[-1,0];d=100;axis([0,100,0,60])fork=1:60Pk=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*')2．动态追击曲线图形程序functiond=chase()Pk=[100,0];P=Pk;Q=[0,0];e=[-1,0];d=100;fork=1:60Pk=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)end3．只有两个点的动态模拟程序functiond=chase()Pk=[100,0];Qk=[0,0];e=[-1,0];fork=1:0.05：60Pk=Pk+0.1*e;Qk=[0,k];e=Qk-Pk;d=norm(e);e=e/d;x=Pk(1);y=Pk(2);u=Qk(1);v=Qk(2);plot(u,v,'o',x,y,'r*',[0,100],[60,0],'og'),pause(.5)end四、有阻力抛射曲线实验对于炮弹飞行，假设阻力与速度成正比，即速度越大则阻力越大。在微分方程中增加阻力项，这是一个比伽力略模型更合符实际的数学模型。)()()()(tykgtytxktx其中，系数k可以通过实验确定。初始条件为：sin)0(,0)0(cos)0(,0)0(00vyyvxx利用初始条件，求解两个二阶常微分方程，得x=dsolve('D2x=-k*Dx','x(0)=0','Dx(0)=198*cos(alfa)')；y=dsolve('D2y=-g-k*Dy','y(0)=0','Dy(0)=198*sin(alfa)')；化简解函数simplify(x)ans=-198*cos(alfa)*(-1+exp(-k*t))/ksimplify(y)ans=-(g*t*k-g-198*sin(alfa)*k+exp(-k*t)*g+198*exp(-k*t)*sin(alfa)*k)/k^2以上结果并不是一个简单的计算公式，计算不简便。考虑用taylor级数展开近似，程序如下symstvgalfakx=dsolve('D2x=-k*Dx','x(0)=0','Dx(0)=v*cos(alfa)');X=taylor(x,3,t)y=dsolve('D2y=-g-k*Dy','y(0)=0','Dy(0)=v*sin(alfa)');Y=simplify(taylor(y,3,t))最后得X=v*cos(alfa)*t-1/2*v*cos(alfa)*k*t^2Y=v*sin(alfa)*t-1/2*g*t^2-1/2*t^2*v*sin(alfa)*k结集号火炮的网上数据为：最大射程2788米。Taylor展开近似，可以较好解释有阻力的抛射曲线模型的合理性，同时也带来计算上的方便。实验程序使用循环语句，试探性地增加飞行时间，以高度数据不大于零为循环终止条件，以确定最大射程，以450做发射角只是一个猜测，最大发射角的研究做为思考题。当阻力系数k=0时就是原来的抛射曲线参数方程，显然不合实际。实验参数（阻力系数）k的取值，类似于猜数游戏。K=0，0.1，0.01，0.02，0.015阻力系数太小，会导致射程太大不合实际；阻力系数太大，导致射程太小不合实验。在一定策略下调整k值，使得射程不断逼近实际数据2788。