常微分方程数值解Word文档下载推荐.docx

1、时间t/s高度h/m速度v/（m/s）加速度a/（m/s2）0.00013.051.0006.56913.1813.292.00026.4226.5613.443.00059.7240.0313.494.000106.553.5013.425.000166.766.8413.256.000240.179.9712.987.000326.592.7712.618.000425.5105.212.159.000536.6117.011.6110.00659.3128.411.0211.00792.9139.110.3812.00937.0149.19.70413.001090.8158.59.01

2、714.001253.4167.28.32515.001424.4175.27.64216.001603.3182.56.98217.001788.9189.26.34718.001980.7195.45.74219.002178.3200.95.17620.002381.4205.94.65921.002589.2210.34.18422.002800.7214.53.73523.003015.6218.33.31824.003234.2221.72.94525.003456.3224.72.62426.003681.9227.22.36727.003909.9229.52.14128.00

3、4139.8231.71.92829.004371.6233.81.73430.004605.3235.71.56731.004841.0237.31.43532.005078.9238.81.34633.005318.21.27634.005558.7241.41.20735.005800.4242.71.14436.006043.3243.91.09037.006287.6245.11.05138.006533.1246.21.03139.006779.7247.21.01940.007027.4248.21.00541.007276.0249.20.99142.007525.7250.3

4、0.97843.007776.4251.30.96944.008028.1252.20.96545.008280.8253.20.96346.008534.5254.20.95847.008789.1255.10.95348.009044.7256.10.94849.009301.2257.00.94350.009558.8258.00.94151.009817.2258.90.93752.0010076.6259.90.93053.0010336.9260.80.92254.0010598.2261.80.91455.0010860.4262.70.90956.0011123.5263.65

5、7.0011387.6264.50.91258.0011652.6265.459.0011918.5266.360.0012185.3267.30.911从表中可以看出，在引擎关闭即燃料燃尽的瞬间，火箭的高度为h=12185.3m，速度v=267.3m/s，加速度a=0.911m/s2。3.3燃料耗尽后物理模型分析当引擎关闭后，火箭继续上升，但由于受重力加速度的影响，速度减小，直至速度减为0时达到最高点。（1）燃料耗尽后火箭质量为320kg。（2）火箭所受合力F=-f-mg;（3）火箭加速度为a=F/m=-f/m-g=;于是可以列出以下微分方程：初值条件为h=12185.3m，v=267.3m

6、/s3.4matlab程序及图像数据记y（1）=h, y（2）=v, y=（y（1）, y（2）T。%-题目3燃料燃尽后函数 rocket M文件-function dx=rocket（t,y）; dx= dy=y（2）;-9.81-0.4*y（2）2/320 %以向量形式表示上述方程组先运行以下程序，以求得v=0时的时间t：%- 题目3燃料燃尽后试验rocket1源程序- %先算出引擎关闭前的速度和高度y0=x（61,1）,x（61,2）; %提取初值for i=1:1500; T=75-0.01*i; %可以先估算出火箭达到最高点的时间小于75s tss=60:0.01:T; t2,y=o

7、de23（rocket,tss,y0）; if y（:,2）=0; %速度小于0时循环终止 break end,1）,b,t2,y（:r）,grid; %绘制全部的高度-时间曲线图4.高度-时间,2）, %绘制全部的速度-时间曲线图5.速度-时间 %计算燃料燃尽后的加速度a2=-9.81-0.4*y（:,2）.2/320;plot（t,a,t2,a2, %绘制全部的加速度-时间曲线图6.加速度-时间v/（m/s2）t2,y（:,1）,y（:,2）,a2;以下为高度、速度、加速度随时间变化的图形：由MATLAB数据计算可知，当t=71.30s 时，火箭上升到最大高度h= 13110m，此时火箭的

8、速度v=0.00309m/s，几乎可认为已经停止，加速度a=-9.81m/s2。3.5 小结全过程中高度、速度、加速度的定性分析：由以上三图的曲线可以看出，起初火箭迅速加速，而且速度值不大时，燃料动力的贡献要超过空气阻力，而质量的减少又会使加速度变大，故而两个因素共同作用使加速度较为稳定。此后速度增加，阻力变大，加速度渐渐变为0。当燃料耗尽时，推动力消失，加速度突变负，而此时速度由增大转为减少。此时速度迅速下降，但是由于它还是正值，故高度上升。速度为0时高度最大，无阻力作用，加速度等于重力加速度。【题目6】一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。已知河水流速v1 与船在

9、静水中的速度v2之比为k。（1） 建立描述小船航线的数学模型，求其解析解；（2） 设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。（3） 若流速v1=0，0.5，1.5，2（m/s），结果将如何？6.1 建立数学模型及分析此种模型的前提是船并不事先知道水速，否则只要调整一个合适的角度，直接沿直线通过即可。现小船不知道水速，则它的策略应为始终使船头瞄准B点。对速度进行xy两个方向的分解，可列出常微分方程如下：其初值条件为x（0）=0,y（0）=-d.6.2 常微分方程解析解两式相除，得变量代换，令, 则有代入有分离变

10、量积分可得代入初值，可解出整理得解析解最终表达式为6.3 MATLAB数值解记x（1）=x, x（2）=y, x=（x（1）, x（2）T。%-题目6函数 guohe.m源文件-function dx=guohe（t,x）v1=1;v2=2;s=（x（1）2+x（2）2）0.5;dx=v1-x（1）*v2/s;-x（2）*v2/s;%-题目6求数值解 guohev.m源文件-x0=0,-100;0.1:70; %经过试验设定终值为70t,x=ode15s（guohe,ts,x0）; %求解方程组plot（t,x）,grid; %画图同时做出x、y方向的位移图1.分位移-时间）x,y/mgtex

11、t（x（t）y（t）plot（x（: %画出航行曲线航行曲线xy,2） %输出数据表得到如下图像： 由图7可以看出，在前30s中，y（t）有较好的的线性，因此，在这个时间段内，小船沿y轴方向的为匀速。从图8可以看出，船合速度的方向为轨迹的切线方向。在刚开始阶段，可以认为传顺水而行，小船在刚开始的速度方向的确变化不大，与上图相吻合。小船在达到x最大值时开始转向，着这一时刻附近的合速度最小。 下面摘取部分数据：t/sx/my/m-100.00000.10000.0999-99.80000.20000.1996-99.60000.30000.2991-99.40000.40000.3984-99.2

12、0000.50000.4975-99.00000.60000.5964-98.80000.70000.6951-98.60000.80000.7936-98.40000.90000.8919-98.20001.00000.9900-98.00001.10001.0879-97.80001.20001.1856-97.60001.30001.2831-97.40001.40001.3804-97.20001.50001.4775-97.00011.60001.5744-96.80011.70001.6710-96.60011.80001.7675-96.400165.10001.5882-0.0

13、98765.20001.4886-0.086665.30001.3889-0.075465.40001.2892-0.064965.50001.1894-0.055265.60001.0896-0.046365.70000.9898-0.038265.80000.8899-0.030865.90000.7900-0.024366.00000.6901-0.018566.10000.5902-0.013566.20000.4902-0.009366.30000.3902-0.005966.40000.2903-0.003366.50000.1903-0.001466.60000.0903-0.0

14、00366.70000.000066.80006.4 解析解与数值解的比较下面是matlab源文件程序：%-题目6 guohebijiao.m源程序-y0=0:-1:-100;x0=50.*（y0./（-100）.0.5-（y0./（-100）.1.5）;,x0,y0,legend（数值解,解析解航行曲线比较由图9可以看出，数值解得到的曲线与解析解的到的曲线几乎完全重合，说明数值解得到的结果还是很准确的。6.5 改变流速（1）v1=0m/s时,更改程序中的数值，可得如下分位移-时间图像及航行曲线：由图像可知，v1=0m/s时，水的流速为0，即小船的合速度等于v2，因此小船沿y轴方向匀速前进。此

15、时小船的位置的解析解为x=0.5*y-0.5*y=0，所以小船过河轨迹为沿y轴方向的一条直线。小船沿y方向的速度，故小船匀速过河。过河所用时间为50s。（2）v1=0.5s时，由数据表可知，t=53.3s时，y=-0.0131,可近似认为已到达岸边。将V1=0.5m/s与v1=1m/s时相比，河水流速减小后，船在x方向上的分速度减小了，因此小船在x方向上行驶的最大值减少。这与途中x最大值减小了一半左右相符。与v1=0时比较，小船过河花的时间并没有增加很多，但路程长了很多，这说明小船的平均合速度大于v2。这是因为实际速度为河水与船水速度的矢量叠加结果。（3）v1=1.5m/s时，由数据表可知，t

16、=114.3s时，x=0.0007m,y=0.000m，可认为小船到达对岸。从图中可以看出，小船在30s，40s时间范围内达到x方向的最大值，然后逆流驶向B点。所以v1越小，船回到B的时间就越短，这一点与常识相符。将解析式求导并令其为零，可以求出x在处，取得最大值，所以k在0,1区间内越大，x取得最大值的点就越靠近x轴。（4）v1=2m/s时利用程序可以求出，小船最终在下游50m处靠岸，不能到达正对岸。也就是说，无论多长时间小船都无法到达正对岸。由解析解也可以证明这一点。，且不恒等于0，所以，对于任意的T0，都有x=，所以小船不可能到达正对岸。将k=1代入得到y=0，所以小船将在y=0时到达x

17、方向的最大值，及当小船到达对岸时，向下游漂的距离最大。由解析式可得，为一条开口向左的抛物线，其轨迹与由matlab作出的轨迹相符。所以小船无法到达B，最终静止在（50,0）处。6.6小结从上述结果可以看出，当静水船速小于等于水速时，不可能通过调整船头始终指向对岸的方法到达对岸。只有当船速大于水速时，才可能用以上方法到达对岸，因为此时静水船速除了可以提供向对岸的分速度外，还能提供逆水流方向的分速度以防船向下漂移。【题目9】两种群相互竞争模型如下：其中x（t），y（t）分别为甲乙两种群的数量；r1,r2为它们的固有增长率；n1,n2为它们的最大容量。s1的含义是，对于供养甲的资源而言，单位数量乙（

18、相对n2）的消耗为单位数量甲（相对n1）消耗的s1倍，对s2可作相应的解释。该模型无解析解，试用数值解法研究以下问题：（1） 设r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值x0=y0=10，计算x（t），y（t），画出它们的图形及相图（x,y）,说明时间t充分大以后x（t），y（t）的变化趋势（人们今天看到的已经是自然界长期演变的结局）。（2） 改变r1,r2,n1,n2,x0,y0,但s1,s2不变（或保持s11），计算并分析所得结果；若s1=1.5（1）,s2=0.7（1）,再分析结果。由此你得到什么结论，请用各参数生态学上的含义作出解释。（3） 试验当s1=0.8（1

19、）, s2=0.7（1）, s2=1.7（1）时又会有什么结果。能解释这些结果吗？9.1模型分析此题的模型很类似教材中的“弱肉强食”模型，只是多考虑了竞争资源的问题。微分方程描述了竞争关系和捕食关系。通过对微分方程组进行数值求解，即得出竞争关系的直观表示。9.2程序代码及输出结果设x（1）=x,x（2）=y。首先编写如下函数：%-题目9函数 jingzh.m源程序-function dx=jingzh（t,x）r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2;dx=r1*x（1）*（1-x（1）/n1-s1*x（2）/n2）;r2*x（2）*（1-s2*x（1）/n1-x

20、（2）/n2）然后运行以下程序%-题目9 jingzhv.m源程序-0.2:15;x0=10,10;t,x=ode23s（jingzh,ts,x0）;t,x;tx,y得到以下图像及数据：txy10.000012.835212.199316.268214.505820.318416.748924.967318.719830.150420.207035.760621.04092.100041.658021.13222.400047.683320.49062.700053.671619.22533.000059.490617.48593.300065.014615.44683.600070.131313.31333.900074.792211.20524.200078.94869.24094.500082.59377.48784.800085.74075.97315.100088.41684.70165.400090.66233.65935.700092.52422.82126.000094.05162.15796.300095.