数学实验——常微分方程数值解.docx

数学实验报告实验4 常微分方程数值解实验4 常微分方程数值解分1 黄浩 2011011743一、 实验目的1. 掌握用MATLAB软件求微分方程初值问题数值解的方法；2. 通过实例学习用微分方程模型解决简化的实际问题；3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。二、 实验内容1. 数学实验第一版（问题2）问题叙述：小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。模型转换及实验过程：（一）从发射到引擎关闭设火箭总质量为m，上升高度为h，瞬时速度为v，瞬时加速度为a，由燃料燃烧时间t=60s，可列如下的方程组：mt=m0-kt=1400-18t ht=vt 其中t0,60 vt=at at=F推 -F阻m(t)-g=32000-0.4v2t1400-18t-9.8因此，上述方程为二元常微分方程组，选择t为自变量，h和v为因变量进行分析。初值条件：h0=0 , v0=0对上述模型，使用ode45()函数求数值解（程序见四.1、四.2），结果如下：xh(t)v(t)a(t)0.000.000.0013.05711.006.5713.1913.30452.0026.4426.5813.45333.0059.7640.0613.49724.00106.5753.5413.43315.00166.7966.8913.26136.00240.2780.0212.98537.00326.7292.8312.61228.00425.79105.2212.15209.00536.99117.1111.616910.00659.80128.4311.021311.00793.63139.1410.380012.00937.85149.189.708313.001091.79158.559.020914.001254.71167.238.330915.001425.93175.227.650216.001604.83182.556.990117.001790.78189.226.359318.001983.13195.275.764619.002181.24200.755.209520.002384.47205.704.694621.002592.36210.184.222022.002804.52214.193.794323.003020.56217.793.412024.003240.08221.013.073025.003462.65223.922.772626.003687.88226.562.504427.003915.58228.972.267728.004145.60231.142.063329.004377.76233.111.889830.004611.86234.911.743331.004847.68236.571.617832.005085.02238.141.506233.005323.85239.611.409534.005564.11240.991.329335.005805.77242.281.265036.006048.72243.501.213937.006292.87244.681.170838.006538.11245.831.130339.006784.48246.961.094740.007031.96248.051.066341.007280.54249.101.045642.007530.19250.121.030843.007780.85251.141.017844.008032.49252.151.002445.008285.12253.160.987646.008538.75254.150.976347.008793.39255.120.969648.009049.01256.070.966349.009305.58257.030.962450.009563.08257.990.952751.009821.52258.950.941252.0010080.93259.900.933753.0010341.30260.830.932854.0010602.62261.750.936355.0010864.86262.670.937056.0011127.98263.610.925857.0011392.04264.540.913858.0011657.03265.460.910659.0011922.96266.350.916160.0012189.78267.260.9170由上表可知，引擎关闭瞬间，火箭的高度为12189.78m，速度为267.26m/s，加速度为0.9170m/s2，火箭至此已飞行60s而高度、速度、加速度随时间的变化曲线如下：(二)从引擎关闭到最高点设引擎关闭时，t1=0，由上一问的结果可知，h0=12189.78m，v0=267.26m/s，m=320kg，则可列二元常微分方程组如下： ht1=vt1 vt1=at1 假设t10,20 at1=-F阻m-g=-0.4v2t1320-9.8因此，可选择t1为自变量，h、v为因变量进行分析（程序见四.3、四.4），实验结果如下：t1ht1vt1at1012190267.26-99.085112416192.70-56.217212585147.43-36.971312716116.10-26.65041282092.80-20.56551290374.30-16.70061296958.95-1479-12.42181306134.05-11.24991309023.21-10.473101310812.99-10.01111131163.11-9.8121213114-6.70-9.8561313102-16.68-10.1481413081-27.15-10.7221513048-38.35-11.6391613004-50.56-12.9951712947-64.36-14.9781812874-80.70-17.9421912784-100.86-22.5162012670-126.76-29.885由上表可知，当t111,12时，vt1有零点，即该区间内某时刻火箭达到最高点。再进行更细致的实验（程序略），设步长为0.01，观察该区间内vt1的零点，如下表所示：t1ht1vt1at111.2613115.350.4175-9.800211.2713115.350.3195-9.800111.2813115.350.2215-9.800111.2913115.360.1235-9.800011.3013115.360.0255-9.800011.3113115.71-0.0411-9.800011.3213115.71-0.1391-9.8000 可以看出，当t1=11.30s，即总时间t=71.30s时，火箭达到最高点，高度为13115.36m，加速度为-9.8m/s2。对0s71s的火箭上升全过程进行作图（程序见四.5）：得出结论：a) 引擎关闭瞬间，火箭的高度为12189.78m，速度为267.26m/s，加速度为0.9170m/s2b) 上升至最高点时，高度为13115.36m，加速度为-9.8m/s2，总时间为71.3s2. 数学实验第一版（问题5）问题叙述：一只小船度过宽为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），结果将如何？模型转换及实验过程： （1）小题. 根据题目要求，小船的指向只能正对B点，也就是说，v2始终由小船指向B点。以A点为原点，水流方向为x轴正方向，AB方向为y轴正方向，建立直角坐标系，如下图所示：根据运动学与动力学公式，不难得到如下的二元常微分方程组：yt=v2cos=v2d-y(d-y)2+x2x(t)=v1-v2sin=v1-v2x(d-y)2+x2v1v2=k为了求出微分方程的解析解，将2式除以1式，得：dxdy=kx2+(d-y)2d-y-xd-y取p=xd-y，则上式可化为：dxdy=k1+p2-p由x=p(d-y)，得dx=ddp-pdy-ydp，代入上式，整理后得：dp1+p2=kdyd-y积分即得：lnp+1+p2=klnCd-y也即：p+1+p2=(Cd-y)k解上面的这个方程,再代入x=0，y=0的值，即得到微分方程的解析解： x=d-y2dk(d-y)k-(d-y)kdk（2）小题.由d=100m，v1=1m/s，v2=2m/s，设计程序，对该情况下微分方程的数值解进行分析（程序见四.6、四.7），结果如下：(省略了前50s的数据)txy5014.4791889.659795113.8356390.801825213.1553591.886015312.4383492.912355411.6859993.880145510.9015994.782805610.0851495.62033579.2384296.39241588.3656297.09313597.4665997.72296606.5450798.28035615.6043298.76148624.6448199.16760633.6722099.49534642.6874999.74418651.6946799.91433660.6972099.9988967-0.0000000181100.00000680.0000000221100.00000690.0000000159100.00000700.0000000724100.00000由上表可见，当出发后67s时，小船已经到达对岸，为了获得更精确的数值解，再设置步长为0.05，细分区间（程序略），结果如下：txy66.500.19720482399.9999208866.550.1472048399.9999576866.600.09720483399.9999843666.650.047204835100.0000066.70-0.0000000114100.0000066.750.0000000395100.0000066.800.0000000073100.0000066.850.0000000739100.0000066.900.0000000697100.00000 由上表可见，当出发66.65s后，小船到达对岸。 然后再给出x-t图像、y-t图像和小船轨迹图：然后对（1）小题获得的解析解，使用matlab作图分析（程序见四.7）：由上图可知，两条轨迹图高度一致，几乎完全重合，这也证明了数值解的正确性。（3）小题.若v1=0m/s ,则显然，小船做匀速直线运动，经过t=d/v2=50s到达终点若v1=0.5m/s ,修改之前的程序（此小题程序见四.8，程序仅给出一种情况，其余略去），得到如下结果：txy485.079692.8849494.381594.4867503.594596.0175512.702197.4529521.679598.7483530.476799.790454-0.00058724100550.0004547610056-0.00019049100通过上面的表格和图片，我们可以得到，当t=54s时，小船到达终点。在进行这一条件下的计算时，matlab的计算已经极为缓慢，结合之前解除的解析解，不难想到这是一个刚性方程，而我一开始使用的是ode45，没有快速处理刚性方程的能力，因此我换用了ode23s。然而在换用函数之后，计算速度仍然十分缓慢，因此，我又使用option函数，将绝对误差从10-6增加到了10-3，这样便可以快速计算。经过比较，降低误差对到达终点的时间的偏离不大。若v1=1.5m/s ，结果如下：txy1102.1716100.00311111.6716100.00081121.171699.99781130.671699.99981140.1716100115-0.000002411001160.000006821001170.00008207100从表格及图像可以看出，当t=115s时，小船到达终点。若v1=2m/s ,则小船的分位移-时间图像和小船轨迹图如下：由上图可见，当v1=2m/s 时，小船无法到达目的地，最终只能到达（50,100）的点附近，即到达距目的地50米的对岸。这是很容易从直观上理解的，因为只有当v2在x轴的分量与v1相抵消，才能使船在x轴方向上偏离后，有能够返回y轴航线的“能力”，显然，此时v2=v1，二者一旦抵消，则v2与v1必等值异号，船便失去了沿y轴前进的动力，即始终停留在原地。当小船行进至（50,100）时便是这种情况，小船始终指向了目标点（0,100），但因v2=v1，丝毫不能前进。通过以上v1=0，0.5，1.5，2（m/s）的四种情况的比较，结合之前v1=1时的讨论，我们可以得出以下两个结论：第一，根据本题的要求，k必须满足：0k1，若k大于1，则因为没有返回航线的能力，始终顺流而下，永远不会到达河的正对岸第二，k越大，到达对岸的时间就越长，而且在沿河方向偏离起始点（或目的地）的最大距离也越大。3. 数学实验第一版（问题9）问题叙述：两种群相互竞争的模型如下：dxdt=r1x(1-xn1-s1yn2)dydt=r2y(1-yn2-s2xn1)其中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）时又会有什么结果。能解释这些结果么？模型转换及实验过程： （1）小题.由于题目中已给出了模型，可以使用matlab直接套用（程序见四.9、四.10，本大题的其余程序均类似，只是修改了参数，在此略去），所得甲乙种群随时间的变化情况以及两种群的相图如下：从以上两图可以看出，当t充分大时，甲种群已成为绝大多数的种群，乙种群已经近乎绝迹。所以随着时间的推移，甲种群的个体数会逐渐增加，直至达到其最大容量100，而乙种群的个体数会逐渐减少，直至灭亡。（2）小题.由于在（1）小题中，乙种群是“弱势”种群，随着时间的推移会逐渐消亡，因此可以通过减小r1、n1、x0或增加r2、n2、y0的手段，来增加乙种群存活的相对能力。我选择了分别增加r2、n2、y0，来依次判断乙种群是否能够得到挽救，结果如下：（为使结构紧凑，省略了相图）调节r2=10时： 由上图可知，当增加r2，即增加种群乙的固有增长率时，虽然种群乙在开始的一段时间占据优势地位，但这并不能挽救种群乙的灭亡。相反，相对于r2=1的情况，当繁殖能力增加10倍之后，种群乙维持的时间更短了。从生态学上分析，在资源一定、没有天敌的情况下，增加繁殖能力对于种群的维持是没有益处的，相反，会因为种群个体数的急剧增加而使个均食物减少，种群内部斗争增加，从而升高死亡率。但是，如果繁殖能力增加了上千倍，那么有可能出现这样的情况，即种群乙在一开始便近乎垂直增长，以至于完全遏制种群甲的繁殖？试调r2=1000：事实是，种群乙绊了自己的脚，存活时间更短了。因为种群乙对甲的竞争能力较弱，因此乙的大量增长对本种群的危害更甚。调节n2=500时： 从上图可以看到，增加n2，即增加种群乙的最大容量，同样无法扭转形势，和未改动之前一样，当t=10左右的时候，种群乙消亡。但是增加最大容量也起到了一点积极作用，即增加了种群乙个体数峰值。调节y0=50时：由上图可见，虽然种群乙仍无法逃脱劫难，但其峰值却大幅提高了，这给我们一点启示，如果继续增加y0，是否会扭转呢？试调y0=200：颓势似乎不可避免，但这时我们很容易能想到，一开始种群乙下降速度太快，是因为超过了最大容量n2=100，种群内竞争激烈，如果再联合调解n2，或许就能及遏制种群甲的繁殖，又不至于对自身造成伤害。联合调节y0=200，n2=300：种群乙终于战胜了种群甲，其原因是在初始时刻，由于大量的乙物种抢食甲的食物，致使甲无法正常生长繁殖，最终会趋于消亡，而因为乙此时还未到最大容量，种内竞争不大，因此仍有正的增长率，这便巩固了种群乙的胜利态势。然后在保持其他参数不变的情况下，调节s1=1.5，s2=0.7，所得结果如下：从上图可见，当乙对甲的竞争力大大增强（甚至强于甲的种内竞争能力）、而甲对乙的竞争力减弱时，乙也能扭转颓势，繁殖到最大容量，并将甲排挤至灭亡的境地。综合上述若干实验，可以得到如下几点结论：1) s为种群的竞争能力，s越大的种群，在与其他种群竞争时，优势越明显，越能排挤其他种群。2) r为种群的固有增长率，即种群的繁殖能力。它只是保证了种群的最高繁殖速度，在没有天敌的情况下，对于提高种群的生存率无益。3) n为种群的最大容量，当种群个体数大于n时，它会抑制种群增加。4) x0（或y0）是初始个体数，初始个体数的差异会造成起始态势的不同，但其之后的演变主要依赖于竞争能力s。当然，对于s较小的种群，在初始值x和最大容量n都很大的时候，也能战胜s较大的种群。（3）小题.当s1=0.8，s2=0.7时，结果如下：因此，在这种情况下，甲乙种群都不会消亡，而是保持相对平衡，而且竞争能力较大的乙的平衡个体数大于甲。下图是甲乙种群的相图：对于结尾一点（45.5,68）附近）的放大图： 从上图看出，甲乙种群最终是处在动态平衡当中，均不会趋于消亡。对于s1=1.5，s2=1.7的情况：由上图可知，在这种情况下，竞争力较差的种群乙渐趋消亡，而种群甲最终占据了主导地位。综合以上两种情况以及（2）小题关于s的讨论，我们可以得出以下结论：1) s代表一种群对另一种群的竞争力。s越大的种群，在与其他种群竞争时，优势越明显，越能排挤其他种群。2) 两种群的s都小于1时，种内竞争为主要因素，因此他们二者不会因为竞争而出现一方消亡、一方完胜的情况，而是最终处于动态平衡，维持在一个相对恒定的数值上。而且s较大的种群，平衡时的个体数也较大。3) 两种群的s都大于1时，种间斗争超过了种内竞争，因此，在其他条件相同的情况下，s更大的一方对另一种群的正常生长繁殖干扰大，随着时间的推移，双方的强弱对比渐趋明显，最终弱势的群体因受到了更大的竞争而出现负增长，最终灭亡。强势的一方也不是一帆风顺，在开始的一个阶段，随着弱势种群的增长，强势方的繁殖速率也会减缓，当弱势种群被排挤出局后，强势一方才恢复了快速的繁殖速率。三、 实验总结本次实验是使用matlab计算常微分方程的数值解，在科研工作中，经常能遇到常微分方程，但是使用解析的方法求出解析解，需要较高的数学技巧，也费时费力，运用matlab求解数值解便是一种快速而准确的方法。通过本次实验的练习，我已经基本掌握了这个方法，能够熟练运用公式，而且能针对实验中出现的问题，进一步展开探究，并予以一一解决。本次三道题都是基于一定的现实背景，不是纯数学意义上的计算，因此，在处理这类问题时，便不能局限于数学的眼光，要灵活运用以前学过的自然科学知识，交叉分析，这样才能使分析结果与实际情况更加吻合，也能够从僵硬的实验数据中获取到更多的有效信息，取得更大的研究成果。四、 程序清单1.第一题第一阶段解常微分方程的函数function dx=ahuojian(t,x)dx=x(2);(32000-0.4*x(2)2)/(1400-18*t)-9.8);end2.第一题从发射到引擎关闭以及绘制图形ts=0:1:60;x0=0,0;t,x=ode45(ahuojian,ts,x0);a=(32000-0.4*x(:,2).2)./(1400-t*18)-9.8);t,x,aplot(t,x(:,1);grid;title(高度-时间曲线);xlabel(t/s);ylabel(h/m);pause;plot(t,x(:,2);grid;title(速度-时间曲线);xlabel(t/s);ylabel(v/m*s-1);pause;plot(t,a);grid;title(加速度-时间曲线);xlabel(t/s);ylabel(a/m*s-2);pause;3.第一题第二阶段解常微分方程的函数function dx=bhuojian(t,x)dx=x(2);(-0.4*x(2)2)/320-9.8;end4.第一题从引擎关闭到最高点ts=0:1:20;x0=12189.78,267.26;t1,x=ode45(bhuojian,ts,x0);a=(-0.4*x(:,2).2)/320-9.8;t1,x,a5.第一题绘制全过程图形ts=0:1:60;x0=0,0;t1,x1=ode45(ahuojian,ts,x0);a1=(32000-0.4*x1(:,2).2)./(1400-t1*18)-9.8);ts=60:1:71;x0=12189.78,267.26;t2,x2=ode45(bhuojian,ts,x0);a2=(-0.4*x2(:,2).2)/320-9.8;plot(t1,x1(:,1),t2,x2(:,1);grid;title(高度-时间曲线);xlabel(t/s);ylabel(h/m);pause;plot(t1,x1(:,2),t2,x2(:,2);grid;title(速度-时间曲线);xlabel(t/s);ylabel(v/m*s-1);pause;plot(t1,a1,t2,a2);grid;title(加速度-时间曲线);xlabel(t/s);ylabel(a/m*s-2);pause;6.第二题小船运行的函数function dx=xiaochuan1(t,x)r=sqrt(100-x(2)2+x(1)2);dx=1-2*x(1)/r;2*(100-x(2)/r;end6.第二题微分方程的数值解ts=0:1:70;x0=0,0;t,x=ode23s(xiaochuan1