《常微分方程matlab求解命方法.ppt》由会员分享，可在线阅读，更多相关《常微分方程matlab求解命方法.ppt（16页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、1/16微分方程与计算机模拟微分方程与计算机模拟 常微分方程数值解方法常微分方程数值解方法捕食者与被捕食者问题捕食者与被捕食者问题有阻力抛射曲线问题有阻力抛射曲线问题卫星轨道模拟问题卫星轨道模拟问题2/16数值方法数值方法求常微分方程初值问题求常微分方程初值问题 求解步骤求解步骤:(1)用用函函数数文文件件定定义义一一阶阶微微分分方方程程(或或方方程程组组)右右端端函函数数；(2)用用MATLAB命令命令ode23()求数值解或绘积分曲线。求数值解或绘积分曲线。使用格式使用格式：T,Y=ode23(F,Tspan,y0)其其中中,Tspan=t0，tN是是常常微微分分方方程程求求解解区区域域，
2、y0是是初初始值始值,F 是包括函数文件名字的符串。是包括函数文件名字的符串。返回值返回值(T,Y)是求解区域内离散数据及对应数值解是求解区域内离散数据及对应数值解。3/16例例7.1.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)命令窗口命令窗口 编辑窗口编辑窗口 4
3、/16常微分方程组初值问题常微分方程组初值问题一阶常微分方程组初值问题数值求解方法一阶常微分方程组初值问题数值求解方法T,y=ode23(F,Tspan,y0)其中其中,F是是函数文件函数文件,表示表示 微分方程右微分方程右端端函数函数Tspan=t0 Tfinal 求解区域求解区域;y0 初始条件初始条件注注:函数函数F(t,y)必须返回列向量必须返回列向量.数值解数值解 y 的每一行对应于列向量的每一行对应于列向量T中的每一行数据中的每一行数据5/16捕食者与被捕食者问题捕食者与被捕食者问题 海岛上有狐狸和野兔海岛上有狐狸和野兔,当野兔数量增多时，狐狸捕食当野兔数量增多时，狐狸捕食野兔导致
4、狐群数量增长野兔导致狐群数量增长;大量兔子被捕食使狐群进入大量兔子被捕食使狐群进入饥饿状态其数量下降饥饿状态其数量下降;狐群数量下降导致兔子被捕食狐群数量下降导致兔子被捕食机会减少机会减少,兔群数量回升。微分方程模型如下兔群数量回升。微分方程模型如下计算计算 x(t)，y(t)当当t0，20时的数据。绘图并分时的数据。绘图并分析捕食者和被捕食者的数量变化规律。析捕食者和被捕食者的数量变化规律。x(0)=100y(0)=20 6/16创建创建MATLAB的函数文件的函数文件function z=fox(t,y)z(1,:)=y(1)-0.015*y(1).*y(2);z(2,:)=-y(2)+0
5、.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/16-兔子数量;-狐狸数量兔-狐数量变化相位图8/16抛射曲线实验抛射曲线实验，假设阻力与速度成正比。在微分方，假设阻力与速度成正比。在微分方程中增加阻力项程中增加阻力项 符号计算方法符号计算方法syms t v g alfa kx=dsolve(D2x=-k*Dx,x(0)=0,Dx(0)=v*cos(alf
6、a);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)9/162008电影集结号展现出视听震撼的战争场面，电影集结号展现出视听震撼的战争场面，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
7、*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;endXmax=xk;plot(x,y,ro)10/16实验数据实验数据:k 0.1 0.01 0.02 0.015Xmax 677.35 3073.15 2433.66 2719.33k=0.02k=0.01511/16嫦娥一号轨道数据实验嫦娥一号轨道数据实验经历四次变轨提速后，卫星才进入地月转移轨道。经历四次变轨提速后，卫星才进入地月转移轨道。第一次变轨卫星由初始轨道进入第一次变轨卫星由初始轨道进入1616小时轨道小时轨道;第二次变轨卫星
8、进入第二次变轨卫星进入2424小时轨道小时轨道;第三次变轨卫星进入第三次变轨卫星进入4848小时轨道小时轨道;第四次变轨卫星进入第四次变轨卫星进入116116小时地月转移轨道。小时地月转移轨道。嫦娥一号卫星进入的初始轨道是周期为嫦娥一号卫星进入的初始轨道是周期为16小时的地球同步轨道。小时的地球同步轨道。卫星进入初始轨道时卫星进入初始轨道时,最大速度大约为最大速度大约为10.3(km/s),而奔月速度需要而奔月速度需要10.9(km/s)。12/16假假设设五五个个轨轨道道上上最最大大速速度度从从10.3(公公里里/秒秒)逐逐步步增增加加到到10.9(公里公里/秒秒)10.3，10.45，10
9、.6，10.75，10.9根据牛顿万有引力定律，地球对卫星的引力大小为根据牛顿万有引力定律，地球对卫星的引力大小为 地球引力参数地球引力参数：GM=3.986005105（km3/s2）卫星运动方程卫星运动方程13/16转换为一阶微分方程组转换为一阶微分方程组初始条件初始条件14/16右端函数的函数文件右端函数的函数文件function z=orbit(t,y)GM=3.986005e05;z(1,:)=y(2);z(2,:)=-GM*y(1)./(y(1).2+y(3).2).(3/2);z(3,:)=y(4);z(4,:)=-GM*y(3)./(y(1).2+y(3).2).(3/2);15/16function Vmax,H=orbitlab(v,h,T)T0=T*60*60;Y0=-(6378+h),v*cos(-pi/2),0,v*sin(-pi/2);T,Y=ode23(orbit,0,T0,Y0);x=Y(:,1);y=Y(:,3);vx=Y(:,2);vy=Y(:,4);V=sqrt(vx.2+vy.2);Vmax=max(V);H=max(x);plot(x,y,0,-(6378+h),0,0,ro)16/16实验数据实验数据(h=200km)Vmax 10.30 10.45 10.60 10.75 10.90 H