[理学]数学实验第二次作业——常微分方程数值求解.doc

实验4常微分方程数值解 实验目的： 练习数值积分的计算； 掌握用MATLAB软件求微分方程初值问题数值解的方法； 通过实例学习用微分方程模型解决简化的实际问题； 了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。 实验内容： 3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。 解答如下： 这是一个典型的牛顿第二定律问题，分析火箭受力情况； 先规定向上受力为正数 建立数学模型： A燃料未燃尽前，在任意时刻（t60s） 火箭受到向上的-F=32000N， 向下的重力G=mg，g=9.8， 向下的阻力f=kv^2, k=0.4, v表示此时火箭速度； 此时火箭收到的合力为F1=（F-mg-f）； 火箭的初始质量为1400kg，燃料燃烧率为-18kg/s; 此刻火箭质量为m=1400-18*t 根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)= (32000-(0.4.*v.^2)-9.8.*(1400-18.*t)) 由此可利用龙格-库塔方法来实现，程序实现如下 Function [dx]=rocket[t,x] %建立名为rocket的方程 m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值 dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)]; %以向量的形式建立方程 [a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式 End; ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点 x0=[0,0]; %输入x的初始值 [t,x]=ode15s(@rocket,ts,x0); %调用ode15s计算 [t,x]; h=x(:,1); v=x(:,2); plot(t,x(:,1)),grid; %绘出火箭高度与时间的关系曲线 title(h-t); xlabel(t/s);ylabel(h/m),pause; plot(t,x(:,2)),grid ; %绘出火箭速度与时间的曲线关系 title(v-t); xlabel(t/s);ylabel(v/m/s),pause; a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t); plot(t,a),grid; %绘出火箭加速度与时间的曲线关系 title(a-t); xlabel(t/s),ylabel(a/m^2/s),pause 火箭高度随时间变化的曲线 火箭速度随时间变化的曲线 火箭加速度随时间变化的曲线 数据过多，故截取部分如下 第一列为时间，第二列为火箭高度，第三列为火箭速度 由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2 B燃料燃尽之后，与A 类似，分析受力如下 火箭受到向上的F=0 向下的重力G=mg，g=9.8， 向下的阻力f=kv^2, k=0.4, v表示此时火箭速度； 此时火箭收到的合力为F2=（-mg-f）； 火箭的初始质量为320kg，恒定 根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320; 程序实现如下 function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数 dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程 ts=60:120; %给出初始时刻，估计终点时刻 x0=[12190,267.26]; %给出x初始值 [t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算 [t,x] plot(t,x(:,1)),grid;