常微分方程的数值解教程文件.doc

常微分方程的数值解精品文档实验4 常微分方程的数值解【实验目的】1掌握用MATLAB软件求微分方程初值问题数值解的方法；2通过实例用微分方程模型解决简化的实际问题；3了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】题3小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。模型及其求解火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式：a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8dh/dt=v又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下：function dx = rocket( t,x )a=(32000-0.4*x(2)2)/(1400-18*t)-9.8;dx=x(2);a;endts=0:1:60;x0=0,0;t,x=ode45(rocket,ts,x0);h=x(:,1);v=x(:,2);a=(32000-0.4*(v.2)./(1400-18*t)-9.8;t,h,v,a;数据如下：thva00013.061.006.5713.1913.302.0026.4426.5813.453.0059.7640.0613.504.00106.5753.5413.435.00166.7966.8913.266.00240.2780.0212.997.00326.7292.8312.618.00425.79105.2212.159.00536.99117.1111.6210.00659.80128.4311.0211.00793.63139.1410.3812.00937.85149.189.7113.001091.79158.559.0214.001254.71167.238.3315.001425.93175.227.6516.001604.83182.556.9917.001790.78189.226.3618.001983.13195.275.7619.002181.24200.755.2120.002384.47205.704.6921.002592.36210.184.2222.002804.52214.193.7923.003020.56217.793.4124.003240.08221.013.0725.003462.65223.922.7726.003687.88226.562.5027.003915.58228.972.2728.004145.60231.142.0629.004377.76233.111.8930.004611.86234.911.7431.004847.68236.571.6232.005085.02238.141.5133.005323.85239.611.4134.005564.11240.991.3335.005805.77242.281.2736.006048.72243.501.2137.006292.87244.681.1738.006538.11245.831.1339.006784.48246.961.0940.007031.96248.051.0741.007280.54249.101.0542.007530.19250.121.0343.007780.85251.141.0244.008032.49252.151.0045.008285.12253.160.9946.008538.75254.150.9847.008793.39255.120.9748.009049.01256.070.9749.009305.58257.030.9650.009563.08257.990.9551.009821.52258.950.9452.0010080.93259.900.9353.0010341.30260.830.9354.0010602.62261.750.9455.0010864.86262.670.9456.0011127.98263.610.9357.0011392.04264.540.9158.0011657.03265.460.9159.0011922.96266.350.9260.0012189.78267.260.92因此，在引擎关闭的瞬间，火箭的速度为267.26m/s，高度为12189.78m，加速度为0.92m/s2。（2）在第二个阶段，火箭的重量保持不变，没有向上的推力，只收到重力和空气阻力。因此有如下关系式：a=dv/dt=(-mg-0.4v2)/m=-0.4v2/320-9.8dh/dt=v假设在80秒之前达到最高点，以60秒时的速度、高度、加速度为初始值进行计算，程序如下：function dx = rocket2( t,x )dx=x(2);-0.4*x(2)2/320-9.8;endts2=60:1:80;x1=12189.78,267.26;t2,x2=ode45(rocket2,ts2,x1);h2=x2(:,1);v2=x2(:,2);a2=-0.4*v2.2./320-9.8;t2,h2,v2,a2;数据如下： t2 h2 v2 a2 60.00 12189.78 267.26 -99.08 61.00 12416.32 192.70 -56.22 62.00 12584.73 147.43 -36.97 63.00 12715.67 116.11 -26.65 64.00 12819.59 92.79 -20.56 65.00 12902.81 74.29 -16.70 66.00 12969.22 58.96 -14.15 67.00 13021.42 45.73 -12.41 68.00 13061.17 33.95 -11.24 69.00 13089.63 23.12 -10.47 70.00 13107.61 12.91 -10.01 71.00 13115.55 3.02 -9.81 72.00 13113.66 -6.80 -9.86 73.00 13101.90 -16.78 -10.15 74.00 13079.96 -27.19 -10.72 75.00 13047.26 -38.34 -11.64 76.00 13002.90 -50.62 -13.00 77.00 12945.47 -64.56 -15.01 78.00 12872.96 -80.96 -17.99 79.00 12782.31 -101.08 -22.57 80.00 12668.88 -127.03 -29.97 可以看到：在第60秒瞬间，加速度发生了突变，从0.92m/s2突变为-99.08m/s2；而在第71秒至第72秒之间，速度从正变为负，即速度反向，说明在第71秒中的某个时刻速度为零，火箭达到了最高点。因此需要对这个时间段进行分析，并且找到速度减小到0的时刻和此时的高度。以0.1为步长，在71s到72s中重新求解微分方程的数值解。71.00 13115.16 2.98 -9.8171.10 13115.41 2.00 -9.8171.20 13115.56 1.02 -9.8071.30 13115.61 0.04 -9.8071.40 13115.57 -0.94 -9.8071.50 13115.42 -1.92 -9.80可见在t=71.3时，速度为0.04，可视为速度为零点，此时最大高度为13115.61，加速度-9.80。综合（1），（2），可以绘出高度，速度和加速度随时间的变化曲线。plot(t,h,t2,h2),xlabel(t/s),ylabel(h/m),title(高度随时间变化曲线);plot(t,v,t2,v2),xlabel(t/s),ylabel(v/(m/s),title(速度随时间变化曲线);aa=a,a2；tt=t,t2；plot(tt,aa),xlabel(t/s),ylabel(a/(m/s2),title(加速度随时间的变化曲线);题6一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比为k。（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；（3）若流速v1=0，0.5，1.5，2m/s，结果将如何。模型及其求解（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B，因此若以B为原点，我们可以得到如下方程组：解初值为（x, y）=( 0, -d) 的常微分方程组，得到解析解为： 其中 c=1/d=0.01,故有事实上，若用matlab中计算微分方程的语句：x,y=dsolve(Dx=v1-v2*x/sqrt(x2+y2),Dy=-v2*y/sqrt(x2+y2),x(0)=0,y(0)=-d,t);则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。（2）d=100m,v1=1m/s,v2=2m/s。求数值解时，由解析解可以看出，此为刚性方程。为避免用ode45s求解时间过长。采用ode23s进行求解。假定100s可以到达对岸。ts=0:0.1:100;x0=0,-100;t,x=ode23s(boat,ts,x0);t,x;plot(t,x);gtext(x);gtext(y);xlabel(时间t/s);ylabel(距离/m);title(x,y与时间t的关系);可以得到数据如下（部分）：t/sx/my/m10.008.93-80.0420.0015.41-60.3630.0018.87-41.4940.0018.71-24.3150.0014.48-10.3450.1014.42-10.2366.500.19-0.0066.600.09-0.0066.700.000可知在t=66.7s时，船到达对岸B点。做 x,y与时间t的关系图：从曲线上可以看出，0到30s这段时间内，y的增长几乎呈线性关系，即小船几乎研直线匀速前进。现在求解析解并将之与数值解对比：function x = jiexi(y,k)x=-0.5*(-0.01)k*y.(k+1)+0.5*(-0.01)(-k).*y.(1-k);endy=-100:1:0;x2=jiexi(y,0.5);plot(x2,y,ro,x(:,1),x(:,2);legend(解析解,数值解,-1);从轨迹曲线中也可以看到，用数值解得到结果与解析解几乎重合，可信度很高。（3）当改变v1的值时，解析式中的k值将发生变化，此时将对结果产生影响。利用MATLAB计算和绘图也可以发现，渡河时间及航行曲线都发生了变化。v1=0时，k=0。说明河水静止不动，这种情况下，小船的总速度就是它在静水中的速度，于是沿着AB直线便可到达对岸，计算结果表示，t=50s时，小船到达B点。v1=0.5m/s时，k=0.25，得到的曲线如下：这种情况下，t=53.33s时，小船到达B点，比起v1=1m/s时，小船在x方向上走过的距离缩短了大约一半，总路程缩短了许多。v1=1.5m/s时，k=0.75，得到的曲线如下：这种情况下， t=114.28s，小船到达对岸，渡河时间明显长了许多。而且由数值解的疏密程度也可以看出，小船花费较少时间久到达x的最大位移处，但是却花了相当长的时间回到x=0的目的地，可见1.5m/s的河水流速给小船到达终点造成了巨大阻碍。v1=2m/s时，k=1，得到的曲线如下：这种情况下，小船无论如何都无法到达B点，只能到达B点下游50米处。从解析式中也可以看到，k=1时，有x(y)=-0.005y2+50。曲线呈开口朝左的抛物线状。从速度合成的三角形上来看，由于v1和v2长度相等，v1的方向也已确定，它们的合速度的方向与v1的夹角不可能大于90，也就是说在x的分位移始终在增加，不可能减小。即使v2沿着水流逆向，也只能使合速度为零，此时正好是小船停在B点下游50米处的情况。类似的问题在高中物理出现过。综合上述曲线，有下图：仔细观察解析式可以发现影响船过河轨迹仅是k值，即船与河水的相对速度，我们不妨作此假设。如果我们用v1=2m/s，v2=4m/s与前面的结果做下对比（我们仅做数值解的对比，因为解析解相同）。结果证明当v1=2m/s，v2=4m/s 时船的航行轨迹与v1=1m/s，v2=2m/s 航行轨迹相同，但时间缩短为33.33s。而且继续实验当v1=4m/s，v2=8m/s时，船的航行轨迹也与前两种情况相同，但过河时间为16.66s，说明过河时间与船速成倒数关系，这是从航行轨迹相同、路程相等可以很自然导出的关系。在实际问题中，人们通常会对水的流速做初步判断，以调整船行驶的方向，而不是单纯的每个时刻都将方向对准目的地。同时由于水的流速也不是始终不变的，在不同的流域和不同的时刻流速都可能不同，本题中只是采取了最为简单的数学模型进行近似计算。题9两种群相互竞争的模型如下：X(t)=r1*x*(1-x/n1-s1*y/n2);Y(t)=r2*y*(1-s2*x/n1-y/n2);其中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，s2=0.7，再分析结果。由此你能得到什么结论，请用各参数生态学上的含义做出解释。（3）实验s1=0.8，s2=0.7和s1=1.5，s2=1.7时会有什么结果，并作出解释。模型及其求解（1）编写程序如下：function dx = compt( 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(2)/n2);endts=0:0.5:20;x0=10,10;t,x=ode23s(compt,ts,x0);t,x;plot(t,x);title(t充分大后x(t),y(t)的变化趋势);xlabel(时间t);ylabel(种群数量);gtext(x(t);gtext(y(t);plot(x(:,1),x(:,2);xlabel(x);ylabel(y);title(相图x,y);设定t从0到20，得出t,x(t),y(t)的数值关系tX(t)Y(t)010.000010.00000.500015.055613.73631.000021.804117.44791.500030.150420.20702.000039.670721.18432.500049.690520.13323.000059.490617.48593.500068.473114.03174.000076.233610.53294.500082.59377.48785.000087.57435.09915.500091.32243.35926.000094.05162.15796.500095.98321.36147.000097.32070.84797.500098.23060.52338.000098.84090.32098.500099.24570.19609.000099.51180.11949.500099.68550.072610.000099.79820.044110.500099.87090.026711.000099.91770.016211.500099.94770.009812.000099.96680.006012.500099.97900.003613.000099.98670.002213.500099.99160.001314.000099.99470.000814.500099.99670.000515.000099.99790.000315.500099.99870.000216.000099.99920.000116.500099.99950.000117.000099.99970.000017.500099.99980.000018.000099.99990.000018.500099.99990.000019.0000100.00000.000019.5000100.00000.000020.0000100.00000.0000由统计数据可以看出，t=17时，乙种群的数量已经为0，之后甲种群的数量达到饱和。由图线更能很好的分析出甲乙两个种群的变化趋势。刚开始时两种群的数量都在增长，但环境的容纳量有限，甲乙两种群不可避免的竞争，于是甲的数量持续增长而乙数量在达到一个峰值后逐渐减少，两种群数量相差悬殊。最后乙灭绝而甲繁荣。究其原因，s1r2，显然甲增加的会更快，结果与（1）相同。而对于r2r1，不妨设r1=1,r2=4,作图。可见开始的时候，凭借着高自然增长率，乙种群数量超过了甲。但之后甲种群数量开始稳步上升，乙种群数量则不可避免的下降，最后甲种群再次繁荣，乙种群灭绝。且达到最后稳定的时间比（1）要短，也许是食物消耗过快的结果。由此可见，乙自然增长率r的提高并不能让乙种群避免灭绝的命运。B. 其次探究最大容量n的影响。显而易见，如果n1n2，则趋势与（1）相同。仅仅是甲的最大数量更大了。对于n2n1，不妨设n2=300,n1=100。作图：可见乙种群最大容量的增大，也不能避免他灭绝的命运。从生态学上来看，这样的结果符合常理。C. 之后探究种群数量初值。对于x0y0，显然趋势和最后结果与（1）相同，甲的增长会更加迅速。对于x0y0，设x0=10,y0=40，作图：可见，乙种群初始数量的增大，虽然在开始时有过一个高峰，但是并不能避免最后灭绝的命运。从生态学上解释，当一个种群生存能力很弱时，即使初始数量很多，最后也会灭绝。下面探究s1，s2。设s1=1.5。s2=0.7。其他因素与（1）相同。可见结果发生了逆转：甲种群最终走向了衰落，几乎灭绝，而乙种群走向了繁荣。由此看出，能否消耗更多供养自己和供养对手的食物和资源，是一个种群能否繁衍生存的关键因素。如果自己的食物和对手的食物都不能有效地“为我所用”，那么一个种群就难逃灭绝的命运。数据如下： tX(t)Y(t)010.000010.00000.500014.164914.87471.000018.809821.18301.500023.178928.67622.000026.402536.81272.500027.943644.98403.000027.766052.68863.500026.226359.66794.000023.815065.84414.500020.981071.24245.000018.062275.92985.500015.262779.98236.000012.702383.45316.500010.447486.41357.00008.505288.91427.50006.868091.00778.00005.508292.74728.50004.392494.18059.00003.486095.35239.50002.755896.303310.00002.171797.070210.50001.707097.685111.00001.339098.175611.50001.048798.565312.00000.820298.873812.50000.640899.117313.00000.500399.309113.50000.390399.459714.00000.304399.577914.50000.237299.670515.00000.184899.742915.50000.144099.799516.00000.112199.843716.50000.087399.878217.00000.068099.905117.50000.052999.926118.00000.041299.942418.50000.032199.955219.00000.025099.965119.50000.019499.972820.00000.015199.9788(3)当s1=0.8，s2=0.7时，两者都小于1。预计到竞争会很激烈，以ts=0:0.5:100;作图：可见，甲乙两个种群竞争十分激烈。开始时增长率几乎相同，到后来乙种群数量超过了甲种群，但是甲种群也并没有因此走向不可逆转的衰亡，最后两者种群数量达到稳定。部分数据如下：tX(t)Y(t)010.000010.00000.500014.771414.86191.000020.769421.07961.500027.529028.26592.000034.231535.65962.500040.081242.46073.000044.590148.11103.500047.705852.46064.000049.650855.64724.500050.736657.92725.000051.251159.56155.500051.392960.74116.000051.313561.61676.500051.121262.30487.000050.852362.84507.500050.562363.30428.000050.260163.69278.500049.959564.03169.000049.671964.34009.500049.397264.617810.000049.137264.866220.000046.365467.380730.000045.677667.986440.000045.506668.136350.000045.466168.171760.000045.457168.179670.000045.455168.181380.000045.454768.181783.000045.454668.181783.500045.454668.181884.000045.454668.1818当t=30以后，两个种群的数量基本保持不变，已经可以视为到达稳定状态。分析：s1=0.8，s2=0.7时，两者都小于1,表示两种群食用供养自己资源的能力要强于食用供养对手资源的能力。在这种情况下，两种群竞争激烈，而没有一方占据绝对优势，所以两种群最后的数量都能够达到稳定值。还可以看出，s微小的差异，会造成最后种群数量较大的不同。若s1=1.5，s2=1.7，两者均大于1。作图如下：很惊奇的，最后甲走向了繁荣而乙走向了灭绝。分析：s1=1.5，s2=1.7时，两者都大于1,表示两种群食用供养对手资源的能力要强于食用供养自己资源的能力。相当于“攻强守弱”。而这种竞争情形下，更强消耗能力成为生存的决定性因素，此种能力更强的一方将会竞争中胜出。与s都小于1的情况截然不同。收集于网络，如有侵权请联系管理员删除