常微分方程数值解.pdf

常微分方程数值解【实验目的】1.练习数值微分的计算；2.掌握用MATLAB软件求微分方程初值问题数值解的方法；3.通过实例学习用微分方程模型解决简化的实际问题；4.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】【题目3小型火箭初始重量为1400kg,其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率 为18kg/s,由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力 正比于速度的平方，比例系数为0.4kg/m,求引擎关闭瞬间火箭的高度、速度、加速度，及 火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。3.1 火箭燃料燃尽前物理模型分析首先设火箭质量为m,高度为h,速度为v,加速度为a,火箭所受合力为F,推力为G,阻力为f,重力加速度为g。(1)由高度、速度、加速度三者之间的关系可得，dh=vdt,dv=adto(2)由火箭燃烧燃料速率为18kg/s可得，火箭质量m(t)=1400；8t(kg)。(3)火箭所受合力 F=G-f-mg=32000-0.4v2-m(t)go(4)火箭加速度为a=m=3鬻=：2-9.8(m/s2)o 1400lot(5)火箭从开始燃烧，到引擎关闭即燃料燃尽所用总时间为t=1080/18=60s。根据上述分析，可以得到一个常微分方程组：(dh=vI dt)dv 32000-0.4v2I -9 81dt 1400-18t其初值条件为：v(0)=0,h(0)=0,t=0;%速度小于。时循环终止 breakend ndplot(t,x(:,l),，r1),grid;%绘制全部的高度-时间曲线 tit 1(，图4.高度-时间，);xlab1(*t/s*);ylab1(1h/m1);pause;plot(t,x(:,2),b ,t2,y(:,2),工)，grid;%绘制全部的速度-时间曲线 title t图5.速度-时间，);xlab1(1t/s,);ylab1(1v/(m/s),);paus;a=(32000-0.4*x(:,2).人2)./(14 0 0-18*七)-9.81;%计算燃料燃尽后的加速度 a2=-9.81-0.4*y(:,2).A2/320;plot(t,a,b ,t2,a2,，r，)，grid;%绘制全部的加速度-时间曲线titlt图6.加速度-时间，);xlab1(*t/s!);ylabel(1v/(m/s2),);t2,y(:,1),y(:,2),a2;%输出数据表以下为高度、速度、加速度随时间变化的图形：图4.高度/寸间300200100图5.速度-时间25010 2030 40 50 60 70 80t/s图6.加速度-时间5020o O-4-6-80-1000 10 20 30 40 50 60 70 80t/s由MATLAB数据计算可知，当t=71.30s时，火箭上升到最大高度h=13U0m,此时火箭的速度v=0.00309m/s,几乎可认为已经停止，加速度a=-9.81m/s2。3.5小结全过程中高度、速度、加速度的定性分析：由以上三图的曲线可以看出，起初火箭迅速加速，而且速度值不大时，燃料动力的贡献 要超过空气阻力，而质量的减少又会使加速度变大，故而两个因素共同作用使加速度较为稳 定。此后速度增加，阻力变大，加速度渐渐变为0。当燃料耗尽时，推动力消失，加速度突 变负，而此时速度由增大转为减少。此时速度迅速下降，但是由于它还是正值，故高度上升。速度为0时高度最大，无阻力作用，加速度等于重力加速度。【题目6一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速V1 与船在静水中的速度V2之比为ko(1)建立描述小船航线的数学模型，求其解析解；(2)设d=100m,Vi=lm/s,v2=2m/s,用数值解法求渡河所需的时间、任意时刻小船的位置 及航行曲线，作图，并与解析解比较。(3)若流速Vi=0,0.5,1.5,2(m/s),结果将如何？6.1 建立数学模型及分析此种模型的前提是船并不事先知道水速，否则只要调整一个合适的角度，直接沿直线通 过即可。现小船不知道水速，则它的策略应为始终使船头瞄准B点。对速度进行xy两个方 向的分解，可列出常微分方程如下：Mx v2xdy v2y其初值条件为x(0)=0,y(0)=-d.6.2 常微分方程解析解两式相除，得色=.虫1三+k 入dy v2y y Jy y变量代换，令u=则有=y+uy dy，dy代入有du i-y=kVu2+l2分离变量积分可得y-k=c(Vu2+I2+u)代入初值，可解出C=d-k整理得解析解最终表达式为X4号尸7)叫6.3 MATLAB数值解用matlab求上述方程组的数值解。记x(l)=x,x(2)=y,x=(x(l),x(2)T0编写M文件如下：%-题目6函数guoh.m源文件-function dx=guoh(t,x)vl=l;v2=2;s=(x(l)A2+x(2)A2)A0.5;dx=vl-x(1)*v2/s;-x(2)*v2/s;nd然后运行以下程序：%-题目6求数值解 guohv.m 源文件-x0=0,-100;ts=0:0.1:70;%经过试验设定终值为70t,x=od15s(guoh，tsf x0);%求解方程组plot(t f x),grid;%画图同时做出x、y方向的位移title(1图1.分位移-时间，)xlab1(*t/s*);ylabel(*x,y/m,);gtxt(*x(t)1);gtext(1y(t)1);paus;plot(x(:,1),x(:,2),grid;%画出航行曲线 title t航行曲线，);xlabel(*x1);ylab1(1y1);t,x(:,l),x(:,2)%输出数据表 得到如下图像:由图7可以看出，在前30s中，y(t)有较好的的线性，因此，在这个时间段内，小船沿y轴 方向的为匀速。从图8可以看出，船合速度的方向为轨迹的切线方向。在刚开始阶段，可以认为传顺水 而行，小船在刚开始的速度方向的确变化不大，与上图相吻合。小船在达到X最大值时开始转向，着这一时刻附近的合速度最小。下面摘取部分数据:t/Sx/my/m00-100.00000.10000.0999-99.80000.20000.1996-99.60000.30000.2991-99.40000.40000.3984-99.20000.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.4001.65.10001.5882-0.098765.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.00146 6.6 0 0 00.0903-0.000366.70000.0000066.80000.000006.4解析解与数值解的比较下面是mat lab源文件程序：%-题目6 guohbi jiao.m源程序-y0=0:-1:-100;x0=50.*(y0./(-100).A0.5-(y0./(-100).A1.5);plot(x(:,l),x(:,2),rlx0,y0,1 b 1)f grid;legend(1数值解I，解析解1)title航行曲线比较D;xlab1(*x1);ylab1(1y1);得到如下图像：0-10-20-30-40-50-60-70-80-90-1000 2 4 6 8 10 12 14 16 18 20 x由图9可以看出，数值解得到的曲线与解析解的到的曲线几乎完全重合，说明数值解得到的结果还是很准确的。6.5 改变流速(1)vl=0m/s时,更改程序中的数值，可得如下分位移-时间图像及航行曲线：分位移出寸间航行曲线U-10IIII II-20-30404-50-60-70-80-90_一-100i-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1x由图像可知，vl=Om/s时，水的流速为0,即小船的合速度等于v2,因此小船沿y轴方 向匀速前进。此时小船的位置的解析解为x=0.5*y-0.5*y=0,所以小船过河轨迹为沿y轴方向的一条直 线。小船沿y方向的速度祟=一品=一詈=V2,故小船匀速过河。过河所用时间为50s.(2)vl=0.5s 时，分位移-时间航行曲线由数据表可知，t=53.3s时，y=0.013L可近似认为已到达岸边。将Vl=0.5m/s与vl=lm/s时相比，河水流速减小后，船在x方向上的分速度减小了，因 此小船在x方向上行驶的最大值减少。这与途中x最大值减小了一半左右相符。与vl=O时比较，小船过河花的时间并没有增加很多，但路程长了很多，这说明小船的平均合速度大于v20这是因为实际速度为河水与船水速度的欠量叠加结果。(3)vl=1.5m/s 时，分位移-时间航行曲线由数据表可知，t=114.3s时，x=0.0007m,y=0.000m,可认为小船到达对岸。从图中可以看出，小船在30s,40s时间范围内达到x方向的最大值，然后逆流驶向B 点。所以vl越小，船回到B的时间就越短，这一点与常识相符。将解析式x=?)_号)叫求导并令其为零，可以求出x在y=-100*2点处,取得最大值，所以k在。1区间内越大，x取得最大值的点就越靠近x轴。(4)vl=2m/s 时分位移-时间航行曲线利用程序可以求出，小船最终在下游50m处靠岸，不能到达正对岸。也就是说，无论 多长时间小船都无法到达正对岸。由解析解也可以证明这一点。=vi-鲁J=V2(1-W身之。，且不恒等于0,dt 1 Jx2+y2 z 7x2+y2/所以，对于任意的T0,都有x=f与也0,所以小船不可能到达正对岸。将匕1代入y=-lOO*?,看得到户0,所以小船将在y=0时到达x方向的最大值，及 当小船到达对岸时，向下游漂的距离最大。由解析式可得,x=捺-0,5|c|y2=50-0.005y2,为一条开口向左的抛物线，其轨迹与由matlab作出的轨迹相符。所以小船无法到达B,最 终静止在(50,0)处。6.6 小结从上述结果可以看出，当静水船速小于等于水速时，不可能通过调整船头始终指向对岸 的方法到达对岸。只有当船速大于水速时，才可能用以上方法到达对岸，因为此时静水船速 除了可以提供向对岸的分速度外，还能提供逆水流方向的分速度以防船向下漂移。【题目9两种群相互竞争模型如下：4)=3(1?戚)(y(t)=r2y(l-s2-)其中x(t),y(t)分别为甲乙两种群的数量；rj2为它们的固有增长率；me为它们的最大 容量。S的含义是，对于供养甲的资源而言，单位数量乙(相对n2)的消耗为单位数量甲(相 对Q)消耗的S倍，对S2可作相应的解释。该模型无解析解，试用数值解法研究以下问题：(1)设1二2=1，ni=n2=100,Si=0.5,s2=2,初值 x=y0=10,计算 x,y(t),画出它们的 图形及相图(x,y),说明时间t充分大以后x(t),y(t)的变化趋势(人们今天看到的已经是自然界 长期演变的结局)。(2)改变1/2川1用2必,丫0,但Si,s2不变(或保持Sil),计算并分析所得结果；若Si=1.5(1)产2=0.7(1),再分析结果。由此你得到什么结论，请用各参数生态学上的含义作出解释。(3)试验当 Si=0.8(1),S2=0.7(1),s2=1.7(l)HtX 会有什么结果。能解释这些结果吗？9模型分析此题的模型很类似教材中的“弱肉强食”模型，只是多考虑了竞争资源的问题。微分方 程描述了竞争关系和捕食关系。通过对微分方程组进行数值求解，即得出竞争关系的直观表 Zj O9.2程序代码及输出结果设x(l)=x,x(2)=y。首先编写如下函数：%-题目9函数jingzh.m源程序-function dx=jingzh(t f x)r1=1;r2=l;nl=100;n2=100;sl=0.5;s2=2;dx=rl*x(1)*(1-x(1)/nl-sl*x(2)/n2);r2*x(2)*(l-s2*x(1)/nl-x(2)/n2)end然后运行以下程序%-题目 9 jingzhv.m 源程序-ts=0:0.2:15;x0=10,10;t,x=oci 23s(jingzhA t s f x0);t,x;plot(t,x),grid;gtext(1x(t)1);gtext(1y(t)1);xlab 1(111);yl北 1(x,y)paus;plot(x(:,l),x(:,2),grid;xlabel(1x1);ylabal(1y1)得到以下图像及数据：tXy0.000010.000010.00000.300012.835212.19930.600016.268214.50580.900020.318416.74891.200024.967318.71981.500030.150420.20701.800035.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.29251.63986.600096.29231.23966.900097.09200.93307.200097.72790.69977.500098.23060.52337.800098.62610.39058.100098.93600.29088.400099.17780.21638.700099.36590.16089.000099.51180.11949.300099.62490.08869.600099.71220.06579.900099.77950.048710.200099.83120.036110.500099.87090.026710.800099.90140.019811.100099.92480.014711.400099.94270.010911.700099.95640.008012.000099.96680.006012.300099.97480.004412.600099.98080.003312.900099.98540.002413.200099.98890.001813.500099.99160.001313.800099.99360.001014.100099.99520.000714.400099.99640.000514.700099.99720.000415.000099.99790.0003X由图像可知，在0,2的时间段内，乙种群数量呈增长状态，之后数量减少，甲种群数量 始终呈上升趋势。在t=9之后，甲乙俩种群数量都趋于稳定。甲种群趋于最大容量(100),乙种群数量趋向0,几乎灭绝。可以认为时间充分长后，甲种群达到最大容量，乙种群几乎灭绝。9.3改变参数L改变rl,r2,即改变固有增长率(1)4 rl=0.3,r2=0.3改变函数中的常数值:%-题目9改变rl r2后的函数-function dx=jingzh(t,x)rl=0.3;r2=0.3;nl=100;n2=100;sl=0.5;s2=2;dx=rl*x(1)*(1-x(1)/nl-sl*x(2)/n2);r2*x(2)*(l-s2*x(1)/nl-x(2)/n2)nd运行以下程序：%-题目9改变rl r2后的程序-ts=0:0.3:30;x0=10,10;t,x=oci23s(jingzh,ts,xO);subplot(2,1,1)r plot(t,x),grid,gtext(1 x(t)1);gtxt(1 y(t)1)r xlabal(111);ylabel(x,y)subplot(2,1,2)rplot(x(:f1),x(:,2),grid;,xlabal(1x1);ylabal(1y1)输出如下图像:由图中数据可以看出，改变二者增长率，最终结果并没有变化。仍然是甲达到100,而 乙彻底灭绝。与之前变化不同的是，达到稳定需要更长的时间，约25年。可以看出，增长 率减少，种群数量变化的速度会减慢。(2)4 rl=5,r2=0.3最终结果仍然没有改变。但不同的是，由于此时甲物种的固有增长率远高于乙，即竞争优势 十分巨大，因此乙物种在初期甚至几乎没有增长就迅速衰减到0.二者约10年时间就达到了 平衡。(3)令门=03/2=5x至此我们基本可以下结论，改变增长率的值不会改变它们最终的演化结果。图中数据表 明，甲物种增长率小，增长速度很慢，相当于减缓了乙物种的灭绝进程。乙初期数量迅速增 长，但是后来随着甲的增多渐渐灭绝。不同的是，乙大约在13年时灭绝，但是甲的数量持 续增长到约23年才稳定。(4)令门=5/2=5150100J 50 0-500 5 10 15 20 25 30可以看出，门、r2的改变不影响最终结果，此时仅用了 2年的时间就达到了最终平衡。而与rl=r2=0.3时的情况相比，x-y图像基本没有改变，说明rl、以同比例改变时，二者之 间的平衡关系不变。2.改变nl,n2,即改变最大容量(1)令nl=500,n2=100,其余数值不变。X由于初期甲物种的数量占最大容量的比例较小(x/nl)即威胁小，所以乙物种增长迅速，最 大值达到了 50左右，但最终仍然灭绝。物种容量的改变并不能影响最终谁会灭绝。x虽然乙的最大容量变大了，但由于甲的作用，乙来不及增长就被灭绝了。乙的最大值与nl=n2=100 时相同。3.改变xO,yO,即改变初值(l).x0=10,y0=100.乙物种的初始数量增大使其灭绝事件稍稍延后，但它灭绝的趋势没有改变。(2).x0=50,y0=10.y(t)x甲的初始数量增大后，乙灭绝的速度更快了。4.小结 从上述的不断控制变量实验中，可以看出，rl和己是甲乙种群的固有增长率，他们决定了 种群在开始一段时间内数量增长的快慢，固有增长率大的种群增长更快。同时决定了到达稳 态的时间快慢，固有增长率大的种群到达稳态时间短。nl和n2是甲乙两种群的最大容量，决定了某一种群最终的稳态值。xO和yO是甲乙两种群的初始值，决定了物种曲线的起点，但不影响最终的稳态。9.4 改变 sls2(1)当 si=0.8(1),s2=0.7(l),s2=1.7(l)时，得到如下图像:小结由上面两个对比试验，可以得出结论。模型中的Is,2s表示的是竞争对手所给予的生存压力。当二者的这一系数都小于1的时候，两个种群会共同生存，在其它条件相似的情况下，系数相对较小的种群会成为数目较大的种群。当其中一个种群的系数大于1而另一种群的系数小于1时，系数小于1的种群会最终 达到最大容量，系数大于1的种群会最终灭亡。当二者的这一系数都大于1的时候，在其它条件相似的情况下，系数相对较小的种群 会最终达到最大容量，系数相对较大的种群会最终灭亡。综上所述si和S2对于两种竞争关系的物种一段时间后数量的变化有着决定作用。si和 S2表示两种物种彼此对资源的竞争关系，物种数量的变化就是由这一因素决定的。竞争力 大的物种最终数量将大于竞争力小的物种，有时两者可以保持一定数量相对稳定较长时间，有时则会出现一个物种数量接近最大容量，另一物种几乎灭绝的情况。所以合理控制si和S2,可以有效保护濒危野生动物，也可以控制一个环境内的生态平 衡。【实验总结】这次实验比第一次要顺利许多。第一次由于对matlab掌握不熟练，经常出现很多自己看不 懂也不知道该怎么改的错误；而经过了上次的多种错误，积累了很多小细节方面的经验，这 次的实验也变得更加顺利。通过本次实验，我掌握用MATLAB求微分方程初值问题数值解的方法，在复习求简单微 分方程解析解的方法的同时，学习了用MATLAB求解微分方程的解析解的方法，并对复杂微 分方程的求解进行了探索。经过这次的实验，不得不感叹数学实验真是奇妙，它大大降低了计算的工作量，为研究 提供了很多的方便。而这些又是建立在建模的基础之上的。