实验4_常微分方程数值解.doc

数学实验 张亚清 2010011805实验4 常微分方程数值解化工系 分0班 2010011805 张亚清【实验目的】1、 掌握用MATLAB软件求微分方程初值问题数值解的方法；2、 通过实例学习用微分方程模型解决简化的实际问题；3、 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】1、 题目3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，为此产生32000N的推力，火箭引擎在燃烧用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。【问题分析及解答】假设火箭上升过程中，重力加速度不随高度而变化，g=9.8m/s2。（1） 关闭引擎前，火箭质量为m，高度为h，速度为v，加速度为a，阻力为f。由题意可知：m=m(t)=1400-18tv=dh/dtf=0.4v2由牛顿第二定律可知： 建立常微分方程组：初值为v=0m/s，h=0m，条件为根据常微分方程组的初值问题，在MATLAB中计算数值解，记x(1)=h，x(2)=v，x=(x(1),x(2)T。具体程序为：function dx=rocket(t,x)dx=x(2);(32000-0.4*x(2)2)/(1400-18*t)-9.8;endts=0:60;x0=0,0;option=odeset(reltol,1e-3,abstol,1e-6);t,x=ode45(rocket,ts,x0,option);t,xplot(t,x(:,1),grid,title(图1.高度-时间)xlabel(t/s)ylabel(h/m)pauseplot(t,x(:,2),grid,title(图2.速度-时间)xlabel(t/s)ylabel(v/(m/s)pausea=(32000-0.4*x(:,2).2)./(1400-18*t)-9.8；plot(t,a),grid,title(图3.加速度-时间)xlabel(t/s)ylabel(a/(m/s2)用MATLAB输出计算得到引擎关闭前的情况：thva0 0.0000 0.0000 13.0571 1 6.5737 13.1887 13.3045 2 26.4444 26.5766 13.4533 3 59.7620 40.0618 13.4972 4 106.5660 53.5352 13.4331 5 166.7909 66.8896 13.2613 6 240.2679 80.0210 12.9853 7 326.7242 92.8288 12.6122 8 425.7867 105.2182 12.1520 9 536.9894 117.1093 11.6169 10 659.7982 128.4334 11.0213 11 793.6301 139.1370 10.3800 12 937.8541 149.1822 9.7083 13 1091.7907 158.5467 9.0209 14 1254.7118 167.2253 8.3309 15 1425.9348 175.2229 7.6502 16 1604.8349 182.5477 6.9901 17 1790.7806 189.2200 6.3593 18 1983.1323 195.2721 5.7646 19 2181.2425 200.7484 5.2095 20 2384.4733 205.7036 4.6946 21 2592.3584 210.1755 4.2220 22 2804.5192 214.1921 3.7943 23 3020.5629 217.7898 3.4120 24 3240.0820 221.0141 3.0730 25 3462.6538 223.9196 2.7726 26 3687.8838 226.5629 2.5044 27 3915.5827 228.9655 2.2677 28 4145.5997 231.1410 2.0633 29 4377.7594 233.1115 1.8898 30 4611.8622 234.9081 1.7433 31 4847.6843 236.5704 1.6178 32 5085.0249 238.1371 1.5062 33 5323.8473 239.6096 1.4095 34 5564.1149 240.9880 1.3293 35 5805.7668 242.2805 1.2650 36 6048.7195 243.5032 1.2139 37 6292.8663 244.6804 1.1708 38 6538.1136 245.8350 1.1303 39 6784.4763 246.9590 1.0947 40 7031.9569 248.0468 1.0663 41 7280.5389 249.0995 1.0456 42 7530.1884 250.1249 1.0308 43 7780.8540 251.1376 1.0178 44 8032.4902 252.1515 1.0024 45 8285.1192 253.1567 0.9876 46 8538.7529 254.1453 0.9763 47 8793.3884 255.1158 0.9696 48 9049.0089 256.0726 0.9663 49 9305.5832 257.0267 0.9624 50 9563.0803 257.9905 0.9527 51 9821.5209 258.9534 0.9412 52 10080.9256 259.9017 0.9337 53 10341.2967 260.8313 0.9328 54 10602.6185 261.7482 0.9363 55 10864.8578 262.6683 0.9370 56 11127.9850 263.6070 0.9258 57 11392.0359 264.5423 0.9138 58 11657.0297 265.4572 0.9106 59 11922.9579 266.3548 0.9161 60 12189.7845 267.2572 0.9170 由此可见，关闭引擎的瞬间，火箭高度h=12190m，速度v=267.26m/s，加速度a=0.9170m/s2。（2） 考虑引擎关闭后的情况，初值为关闭引擎瞬间的数值，条件是t60s，且速度v0，用matlab找出符合上述条件的时间范围60，T。记y(1)=h，y(2)=v，y=(y(1),y(2)T程序如下：function dy=rocket2(t,y)dy=y(2);-0.4*y(2).2/320-9.8;endformat short gts=0:60;x0=0,0;option=odeset(reltol,1e-3,abstol,1e-6);t,x=ode45(rocket,ts,x0,option); for n=1:50 T=71.5-0.01*n; tss=60:0.02:T; y0=x(61,1),x(61,2);%获得初值 option=odeset(reltol,1e-3,abstol,1e-6); t2,y=ode45(rocket2,tss,y0,option); t2,y; if y(:,2)=0 break endendplot(t,x(:,1),b,t2,y(:,1),r),grid,title(图1.高度-时间)xlabel(t/s)ylabel(h/m)pauseplot(t,x(:,2),b,t2,y(:,2),r),grid,title(图2.速度-时间)xlabel(t/s)ylabel(v/(m/s)pausea=(32000-0.4*x(:,2).2)./(1400-18*t)-9.8;a2=-9.8-0.4*y(:,2).2/320;plot(t,a,b,t2,a2,r),grid,title(图3.加速度-时间)xlabel(t/s)ylabel(a/(m/s2)A=t2,y,a2途中蓝色曲线表示引擎关闭前状态，红色表示引擎关闭后状态。由MATLAB输出表格，关闭引擎后的状态为：thva60.0012189.784523267.25716264-99.0829887360.0212195.109937265.28865212-97.77258617760.0412200.396242263.34608064-96.48894773560.0612205.64395261.42892044-95.23135055360.0812210.853565259.53665778-93.99909591460.1012216.025579257.66879237-92.79150820460.1212221.160476255.82483739-91.60793428360.1412226.258748254.00429775-90.44772909460.1612231.320868252.20670292-89.310276245.71.1413115.1313591.5880446369-9.803152357271.1613115.1611591.3919890565-9.802422041971.1813115.1870381.1959471165-9.801787861971.2013115.2089970.99991689453-9.801249792271.2213115.2270350.80389646942-9.800807811971.2413115.2411530.6078839208-9.800461903671.2613115.251350.41187732928-9.800212053771.2813115.2576280.21587477646-9.800058252471.3013115.2599850.019874344877-9.8000004937当t=71.30时，火箭上升到最大高度h=13115m，此时火箭的速度v=0.019874m/s，近似为0，加速度a=-9.8m/s2。【结果分析与讨论】由图像可以明显看出，关闭引擎前，火箭高度和速度均随时间增长而增长，加速度在开始有小幅上升，达到最大值后开始下降，下降速度先快后慢。关闭引擎后，高度仍在一定时间范围内持续增长，速度大幅减小直到为零，加速度瞬间降为-99m/s2，之后逐渐增加，趋于-g。2、 题目6一只小船渡过宽为d的河流，目标是起点A正对着的B点，已知河水流速v1与船在静止的水中的速度v2之比为k。（1） 建立描述小船航线的微分方程模型，求其解析解；（2） 设d=100m，v1=1m/s，v2=2m/s，用数值方法求渡河所需时间、任意时刻小船位置及航行曲线，作图，并与解析解比较；（3） 若流速v1=0,0.5,1.5,2（m/s），结果将如何。【问题分析】 如果小船驾驶者知道水流的速度，则，只要使v10，都有，所以小船不可能到达B点。将k=1代入得到y=0，所以小船将在y=0时到达x方向的最大值，及当小船到达对岸时，向下游漂的距离最大。由解析式得出，为一开口向左的抛物线，其轨迹与由MATLAB做出的轨迹相符。所以小船无法到达B点，最终静止在（50,0）处。3、 题目9两种群相互竞争模型如下：其中x(t)，y(t)分别为甲乙两种群的数量；r1，r2为它们的固有增长率；n1，n2为它们的最大容量。s1的含义是，对于供养甲的资源而言，单位数量乙（相对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)，s2=0.7(1)，s2=1.7(1)时又会有什么结果。能解释这些结果吗？【问题分析及解答】r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值x0=y0=10。由生态学上分析甲和乙两种群是竞争关系。为考察模型中甲乙两种群密度随时间的演变过程，用MATLAB求方程的数值解，记x(1)=x，x(2)=y，x=(x(1),x(2)T。程序如下：function xdot=shier(t,x)r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2;xdot=diag(r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-s2*x(1)/n1-x(2)/n2)*x;endts=0:0.1:15;x0=10,10;t,x=ode45(shier,ts,x0);t,xplot(t,x),grid,gtext(x(t),gtext(y(t),pause,plot(x(:,1),x(:,2),grid,xlabel(x),ylabel(y)得到两种群密度随时间的变化： tx(t)y(t)010.0000 10.0000 0.110.8805 10.7120 0.211.8235 11.4454 0.312.8309 12.1962 0.413.9044 12.9594 0.515.0453 13.7294 0.616.2544 14.4999 0.717.5321 15.2651 0.818.8783 16.0173 0.920.2924 16.7470 121.7729 17.4454 1.123.3182 18.1043 1.224.9260 18.7158 1.326.5933 19.2731 1.428.3169 19.7697 1.530.0927 20.1998 1.631.9163 20.5581 1.733.7827 20.8401 1.835.6862 21.0418 1.937.6212 21.1597 239.5833 21.1919 2.141.5662 21.1405 2.243.5633 21.0083 2.345.5684 20.7987 2.447.5754 20.5155 2.549.5788 20.1628 2.651.5729 19.7455 2.753.5525 19.2687 2.855.5126 18.7382 2.957.4486 18.1600 359.3560 17.5410 3.161.2306 16.8880 3.263.0684 16.2089 3.364.8658 15.5115 3.466.6195 14.8040 3.568.3272 14.0911 3.669.9862 13.3773 3.771.5939 12.6673 3.873.1483 11.9655 3.974.6479 11.2756 476.0915 10.6011 4.177.4785 9.9450 4.278.8087 9.3099 4.380.0824 8.6981 4.481.3002 8.1114 4.582.4629 7.5510 4.683.5703 7.0171 4.784.6233 6.5102 4.885.6228 6.0303 4.986.5701 5.5775 587.4663 5.1514 5.188.3129 4.7517 5.289.1114 4.3777 5.389.8636 4.0284 5.490.5712 3.7030 5.591.2363 3.4001 5.691.8606 3.1186 5.792.4458 2.8577 5.892.9938 2.6162 5.993.5063 2.3932 693.9851 2.1877 6.194.4320 1.9985 6.294.8486 1.8245 6.395.2367 1.6647 6.495.5979 1.5178 6.595.9338 1.3829 6.696.2460 1.2592 6.796.5358 1.1461 6.896.8047 1.0427 6.997.0540 0.9483 797.2849 0.8623 7.197.4987 0.7839 7.297.6966 0.7123 7.397.8795 0.6470 7.498.0487 0.5874 7.598.2050 0.5331 7.698.3493 0.4837 7.798.4824 0.4388 7.898.6052 0.3981 7.998.7184 0.3611 898.8228 0.3275 8.198.9189 0.2970 8.299.0074 0.2692 8.399.0890 0.2439 8.499.1641 0.2210 8.599.2331 0.2003 8.699.2965 0.1814 8.799.3549 0.1644 8.899.4085 0.1489 8.999.4577 0.1349 999.5030 0.1222 9.199.5446 0.1106 9.299.5828 0.1002 9.399.6178 0.0907 9.499.6500 0.0821 9.599.6795 0.0744 9.699.7065 0.0673 9.799.7313 0.0610 9.899.7541 0.0552 9.999.7749 0.0500 1099.7941 0.0452 10.199.8116 0.0409 10.299.8276 0.0370 10.399.8423 0.0335 10.499.8558 0.0304 10.599.8681 0.0275 10.699.8794 0.0249 10.799.8898 0.0225 10.899.8992 0.0204 10.999.9079 0.0184 1199.9158 0.0167 11.199.9231 0.0151 11.299.9297 0.0137 11.399.9358 0.0124 11.499.9413 0.0112 11.599.9464 0.0101 11.699.9510 0.0092 11.799.9553 0.0083 11.899.9592 0.0075 11.999.9627 0.0068 1299.9659 0.0061 12.199.9689 0.0056 12.299.9716 0.0050 12.399.9741 0.0046 12.499.9763 0.0041 12.599.9784 0.0037 12.699.9803 0.0034 12.799.9820 0.0031 12.899.9836 0.0028 12.999.9850 0.0025 1399.9863 0.0023 13.199.9875 0.0020 13.299.9886 0.0019 13.399.9896 0.0017 13.499.9905 0.0015 13.599.9914 0.0014 13.699.9921 0.0012 13.799.9928 0.0011 13.899.9935 0.0010 13.999.9940 0.0009 1499.9946 0.0008 14.199.9950 0.0008 14.299.9955 0.0007 14.399.9959 0.0006 14.499.9962 0.0006 14.599.9966 0.0005 14.699.9969 0.0005 14.799.9972 0.0004 14.899.9974 0.0004 14.999.9976 0.0003 1599.9978 0.0003 由图像可知，在0,2的时间段内，乙种群数量呈增长状态，之后数量减少，甲种群数量始终呈上升趋势。在t=9之后，甲乙俩种群数量都趋于稳定。甲种群趋于最大容量（100），乙种群数量趋向0，几乎灭绝。可以认为时间充分长后，甲种群达到最大容量，乙种群几乎灭绝。改变r1,r2,n1,n2,x0,y0，但s1,s2不变。当r1=2，其他参数不变，得到如下图像：由此可见，甲乙两种群的数量变化情况大致不变，只是比之前更快到达稳态，甲种群比乙种群先达到稳态，即甲达到最大容量，乙几乎灭绝。当r2=2，其他参数和（1）中一样，得到如下图像：可以看出，总体数量变化趋势和（1）中情况相同，乙比甲先达到稳态，在前期乙的增长更快。比较分别改变r1和r2的情况，固有增长率大的种群在开始阶段数量增长更快，也更早到达稳定状态。当n1=200，其他条件同（1）时，得到如下图像：由图像可知，甲种群趋向于最大容量（200），乙种群仍是趋向于0，乙种群开始增长速度小于甲种群，两者几乎同时到达稳态。当n2=200，其他条件同（1）时，得到如下图像：可以看到，甲乙的稳态和（1）一样，乙仍是几乎灭绝。由此可见，甲种群数量趋向于最大容量（n1），乙种群数量趋向于0。当x0y0，如x0=20，y0=10时，其他参数同（1）中，得到如下图像：当x01)，s2=0.7(1)时，其他数据同（1）中，得到如下图像：由图中可见，长时间后乙种群趋向最大容量，甲种群几乎灭绝。【结果分析与讨论】从上述的不断控制变量实验中，可以看出，r1和