数学实验报告-利用MALTAB计算常微分方程数值解.docx

《数学实验报告-利用MALTAB计算常微分方程数值解.docx》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、实验二常微分方程数值解一、火箭飞行器㈠问题描述小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭，设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度及火箭达到最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。㈡方法与公式1、简要分析本题的求解需要用到常微分方程，而整个过程又被分为两个阶段：火箭加速上升阶段和燃料燃尽后减速的阶段。由题目易知第一个阶段持续时间T1=108018=60s列出第一阶段的方程组：设M0为火箭本身质量

2、，m为燃料质量，g为重力加速度=9.8m/s2,燃料燃烧率为a，空气阻力的比例系数为k，F为推进力。M0=1400-1080=320kg;v=F-kv2-M0+mgM0+mm=-a初值v=0，m=1080。由以上各式可以求出t=T1时火箭的速度。再求解第二阶段：v=-kv2-M0gM0m=0可以求出火箭速度降为0的时刻。将整个过程中的时间向量以及速度向量联合起来，利用第三章所学插值与数值积分的方法可以求得任意时刻火箭的近似高度。2、方法求解常微分方程时，我分别采用了自己编写的欧拉公式、改进欧拉公式、4级4阶龙格-库塔公式，以及MATLAB自带的龙格-库塔方法，求解数值积分时采用辛普森公式。由

3、于Matlab自带的Simpson公式是自适应的，因此需要使用自己在上一次实验时所编的Simpson公式。㈢结果与分析1、各种公式的对比首先，我作出了各种不同公式计算得到的火箭速度随时间变化的图像，图如下：从图中可以看出，各种公式计算得到的结果基本一致，为确定其区别，将图像放大，放大约2000倍后，得到下图：分析：从图中可以看出，自编欧拉公式距离MATLAB自带龙格-库塔公式最远，精度最差；自编的改进欧拉公式和自编的龙格-库塔公式结果基本一致，两者中自编龙格-库塔公式距MATLAB自带龙格-库塔公式的结果稍近。与之前的分析基本一致。然而产生自编龙格-库塔公式与MATLAB自带龙格-库塔公式之

4、间的差距的原因还未知。由于MATLAB自带龙格-库塔公式精度较高，因此以下各项实验不再计算其它几项公式的结果。2、第一阶段火箭关闭瞬间的速度v=267.261240773261m/s关闭前瞬间的加速度a=F-M0g-kv2M0=0.914286475421320m/s2此时火箭的高度h=12189.7791507305m3、第二阶段由初始速度以及常微分方程可以求得火箭达到最高点的时间约为71.3s；此时火箭的高度h=13115.7148739887m加速度a=-9.80000406052111m/s24、整体过程下图为火箭加速度与时间的关系。由图可以看出，火箭一开始的加速度很大，随着时间的推

5、移，火箭的燃料有所减少，与此同时速度有所上升。两者中前者使火箭的加速度增大，后者使其减小，综合作用，最终体现为加速度先略有上升，然后慢慢减小。当火箭中燃料燃尽时，火箭丧失推动力，因而加速度急剧减小为负值。此后火箭速度不断减小，导致火箭所受阻力逐渐减小，因而加速度的绝对值有所减小，直到最终火箭速度降为零，火箭不受阻力，仅受重力，加速度为重力加速度。下图为火箭速度与时间的关系：此图可以看作是由第一幅图对时间积分所得结果，本图的斜率对应第一幅图的值。第一阶段火箭加速度为正，因此速度不断增加，只是增加的速度不断减慢。燃料燃尽后，加速度变为负值，因此速度开始急剧下降，与此同时下降的速率不断减小。最终火

6、箭速度降为0。下图为火箭高度与实践的关系：此图可以看作是由第二幅图对时间积分所得结果，本图的斜率对应第一幅图的值。一开始或减速度较小，因此高度缓慢增加，之后增加的速度不断提升，直到火箭的燃料耗尽，此后速度不断减小，但仍为正值，因此火箭继续向上飞行，只是高度提升的速度逐渐变慢。知道最终火箭速度降为零。3、总结这是自己接触数学实验以来第一次解决相对复杂的实际问题，总体上还是有些不适应。虽然最终得以成功解决，但是编程的逻辑以及脚本之间的层次没有考虑清楚，需要在之后的学习中不断加强相关能力。当然也有不小收获：(1)了解了MATLAB解决常微分方程组的原理与步骤；(2)进一步熟悉了MATLAB的相关操

7、作，包括画图、调试、编写脚本等；(3)加强了对之前插值与数值积分的理解，学习将其应用到实际问题的解决中。㈣程序清单1、第一阶段的微分方程组functiondx=fly(t,y)M0=320;a=18;k=0.4;g=9.8;F=32000;dx=[(-k*y(1)^2+F-(M0+y(2))*g)/(M0+y(2));-a*sign(y(2))];2、第二阶段微分方程组functiondx=fly1(t,y)