实验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];数据如下：t h v a0 0 0 13.061.00 6.57 13.19 13.302.00 26.44 26.58 13.453.00 59.76 40.06 13.504.00 106.57 53.54 13.435.00 166.79 66.89 13.266.00 240.27 80.02 12.997.00 326.72 92.83 12.618.00 425.79 105.22 12.159.00 536.99 117.11 11.6210.00 659.80 128.43 11.0211.00 793.63 139.14 10.3812.00 937.85 149.18 9.7113.00 1091.79 158.55 9.0214.00 1254.71 167.23 8.3315.00 1425.93 175.22 7.6516.00 1604.83 182.55 6.9917.00 1790.78 189.22 6.3618.00 1983.13 195.27 5.7619.00 2181.24 200.75 5.2120.00 2384.47 205.70 4.6921.00 2592.36 210.18 4.2222.00 2804.52 214.19 3.7923.00 3020.56 217.79 3.4124.00 3240.08 221.01 3.0725.00 3462.65 223.92 2.7726.00 3687.88 226.56 2.5027.00 3915.58 228.97 2.2728.00 4145.60 231.14 2.0629.00 4377.76 233.11 1.8930.00 4611.86 234.91 1.7431.00 4847.68 236.57 1.6232.00 5085.02 238.14 1.5133.00 5323.85 239.61 1.4134.00 5564.11 240.99 1.3335.00 5805.77 242.28 1.2736.00 6048.72 243.50 1.2137.00 6292.87 244.68 1.1738.00 6538.11 245.83 1.1339.00 6784.48 246.96 1.0940.00 7031.96 248.05 1.0741.00 7280.54 249.10 1.0542.00 7530.19 250.12 1.0343.00 7780.85 251.14 1.0244.00 8032.49 252.15 1.0045.00 8285.12 253.16 0.9946.00 8538.75 254.15 0.9847.00 8793.39 255.12 0.9748.00 9049.01 256.07 0.9749.00 9305.58 257.03 0.9650.00 9563.08 257.99 0.9551.00 9821.52 258.95 0.9452.00 10080.93 259.90 0.9353.00 10341.30 260.83 0.9354.00 10602.62 261.75 0.9455.00 10864.86 262.67 0.9456.00 11127.98 263.61 0.9357.00 11392.04 264.54 0.9158.00 11657.03 265.46 0.9159.00 11922.96 266.35 0.9260.00 12189.78 267.26 0.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 a260.00 12189.78 267.26 -99.0861.00 12416.32 192.70 -56.2262.00 12584.73 147.43 -36.9763.00 12715.67 116.11 -26.6564.00 12819.59 92.79 -20.5665.00 12902.81 74.29 -16.7066.00 12969.22 58.96 -14.1567.00 13021.42 45.73 -12.4168.00 13061.17 33.95 -11.2469.00 13089.63 23.12 -10.4770.00 13107.61 12.91 -10.0171.00 13115.55 3.02 -9.8172.00 13113.66 -6.80 -9.8673.00 13101.90 -16.78 -10.1574.00 13079.96 -27.19 -10.7275.00 13047.26 -38.34 -11.6476.00 13002.90 -50.62 -13.0077.00 12945.47 -64.56 -15.0178.00 12872.96 -80.96 -17.9979.00 12782.31 -101.08 -22.5780.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 点。
已知河水流速v 1与船在静水中的速度v 2之比为k 。
（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m ，v 1=1m/s ，v 2=2m/s ，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；（3）若流速v 1=0，0.5，1.5，2m/s ，结果将如何。
模型及其求解（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B ，因此若以B 为原点，我们可以得到如下方程组：解初值为（x, y ）=( 0, -d) 的常微分方程组，得到解析解为：其中 c=–1/d=–0.01,故有1v x dxv dt v y dy dt ⎧=-⎪⎪⎨⎪=-⎪⎩1111022k k kkx c ycy+--+-=1111()(0.01)(0.01)22kk kkx y yy+--=--+-事实上，若用matlab中计算微分方程的语句：[x,y]=dsolve('Dx=v1-v2*x/sqrt(x^2+y^2)','Dy=-v2*y/sqrt(x^2+y^2)','x(0) =0','y(0)=-d','t');则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。