54实验四_常微分方程数值解

实验四常微分方程数值解材料系材932009011976邓陟一．实验目的1.用MATLAB软件掌握求微分方程数值解的方法，并对结果作初步分析。2.通过实例学习用微分方程模型解决简化的实际问题。二．实验内容1.小型火箭初始质量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭达到最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。可以根据牛顿第二定律Fma列出火箭不同飞行时期满足的运动方程。1)燃料未用尽前(060t)火箭飞行时受到燃料燃烧产生的推力、自身的重力以及空气阻力（向下）。列出微分方程200=00,00pdhvdtFdvkvgadtmutmuthv其中pF为燃料燃烧产生的推力，g为重力加速度（取9.8m/s2），k为空气阻力比例系数，0m为火箭初始质量，u为燃料燃烧率。编写M文件rocket1.mfunctiondx=rocket1(t,x)Fp=32000;g=9.8;k=0.4;m0=1400;u=18;dx=[x(2);Fp/(m0-u*t)-g-k*x(2)^2/(m0-u*t)];编写程序ts=0:0.1:60;x0=[0,0];[t,x]=ode45(@rocket1,ts,x0)得到引擎关闭瞬间火箭的高度6012190hm，速度60267.26vm/s，带入微分方程算得加速度0.92am/s2。2)燃料用尽瞬间至火箭达到最高点火箭飞行时受到自身的重力以及空气阻力（向下）。列出微分方程2=6012190,60267.26dhvdtdvkvgadtmhv其中g为重力加速度（取9.8m/s2），k为空气阻力比例系数，m为燃料用尽后的火箭质量。编写M文件rocket2.mfunctiondx=rocket2(t,x)g=9.8;k=0.，4;m=320;dx=[x(2);-g-k*x(2)^2/m];编写程序ts=60:0.1:75;x0=[12190,267.26];[t,x]=ode45(@rocket2,ts,x0)算得在71.3ts时0v，此时火箭的高度71.313116hm，加速度9.8agm/s2。3)最高点变加速下落到最终匀速下落至地面火箭飞行时受到自身的重力以及空气阻力（向上）。列出微分方程2=71.313116,71.30dhvdtdvkvgadtmhv其中g为重力加速度（取9.8m/s2），k为空气阻力比例系数，m为燃料用尽后的火箭质量。为画出火箭高度、速度、加速度随时间变化的图形，须描述火箭整个的运动过程。编写M文件rocket.mfunctiondx=rocket(t,x)Fp=32000;g=9.8;k=0.4;m0=1400;u=18;m=320;ift=60dx=[x(2);Fp/(m0-u*t)-g-k*x(2)^2/(m0-u*t)];endift60&t=71.3dx=[x(2);-g-k*x(2)^2/m];endift71.3dx=[x(2);-g+k*x(2)^2/m];end编写程序ts=0:0.1:225.2;x0=[0,0];[t,x]=ode45(@rocket,ts,x0);fori=1:2253ift(i)=60a(i)=(32000-0.4*x(i,2).^2)/(1400-18*t(i))-9.8;endift(i)60&t(i)=71.3a(i)=-0.4*x(i,2).^2/320-9.8;endift(i)71.3a(i)=0.4*x(i,2).^2/320-9.8;endendsubplot(3,1,1),plot(t,x(:,1));subplot(3,1,2),plot(t,x(:,2));subplot(3,1,3),plot(t,a);高度、速度、加速度随时间变化的图形2.一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速1v与船在静水中的速度2v之比为k。(1)建立描述小船航线的数学模型，求其解析解；(2)设d100m，1v1m/s，2v2m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行路线，作图，并与解析解比较；(3)若流速1v0,0.5,1.5,2(m/s)，结果将如何。1)建立数学模型以B点为坐标系原点建立直角坐标系，A点坐标为0,d。设小船航线上任意一点的坐标为,xy，2v与1v夹角为，列出微分方程1212222222cossin(0)0,0xydxxvvvvvdtxydyyvvvdtxyxyd两式相除，得22xydxxkdyyy令xpyy，上式可化简为21dpykpdy根据边界条件0pd积分得到21kyppd12kkyypdd112kkdyyxdd此即小船航线的解析表达式。2)用数值法求解编写M文件boat.mfunctiondx=boat(t,x)v1=1;v2=2;dx=[v1-v2*x(1)/sqrt(x(1)^2+x(2)^2);-v2*x(2)/sqrt(x(1)^2+x(2)^2)];编写程序ts=0:70;x0=[0,-100];[t,x]=ode45(@boat,ts,x0);yt=-100:0;xt=50*((yt/(-100)).^0.5-(yt/(-100)).^1.5);plot(x(1:68,1),x(1:68,2),'r',xt,yt,'b')当67ts时，61.11100x，61.33100y，可以认为此时小船已到达B点，则渡河时间67ts。任意时间小船位置txytxy00.0000-100.00003419.2383-34.361910.9899-98.00003519.2418-32.627921.9595-96.00033619.2062-30.916932.9082-94.00093719.1309-29.230743.8356-92.00233819.0155-27.570654.7414-90.00453918.8596-25.938565.6251-88.00794018.6632-24.335976.4861-86.01274118.4256-22.764787.3241-84.01934218.1468-21.227198.1386-82.02804317.8264-19.7249108.9290-80.03904417.4646-18.2603119.6949-78.05274517.0614-16.83531210.4358-76.06964616.6171-15.45181311.1511-74.09004716.1325-14.11151411.8403-72.11434815.6076-12.81711512.5028-70.14304915.0425-11.57181613.1382-68.17655014.4376-10.37831713.7457-66.21545113.7939-9.23901814.3249-64.26025213.1128-8.15561914.8752-62.31155312.3962-7.12942015.3960-60.36975411.6455-6.16242115.8866-58.43575510.8615-5.25842216.3464-56.51015610.0457-4.42012316.7748-54.5935579.2004-3.64922417.1713-52.6868588.3280-2.94642517.5351-50.7907597.4303-2.31512617.8657-48.9061606.5097-1.75752718.1624-47.0338615.5692-1.27362818.4246-45.1749624.6113-0.86622918.6516-43.3303633.6389-0.53583018.8429-41.5012642.6550-0.28383118.9977-39.6886651.6626-0.11093219.1156-37.8937660.6649-0.01773319.1960-36.1177670.00000.0000小船航行路线（数值解与解析解）可以看出数值解与解析解重合得很好。3)对于不同流速①10v修改M文件boat.mfunctiondx=boat(t,x)v1=0;v2=2;dx=[v1-v2*x(1)/sqrt(x(1)^2+x(2)^2);-v2*x(2)/sqrt(x(1)^2+x(2)^2)];编写程序ts=0:50;x0=[0,-100];[t,x]=ode45(@boat,ts,x0);yt=-100:0;xt=0*yt;plot(xt,yt,'b',x(:,1),x(:,2),'r')水速为0时船作匀速直线运动，航行路线为直线，50ts时到达对岸B点。任意时刻小船位置txytxy00-100250-5010-98260-4820-96270-4630-94280-4440-92290-4250-90300-4060-88310-3870-86320-3680-84330-3490-82340-32100-80350-30110-78360-28120-76370-26130-74380-24140-72390-22150-70400-20160-68410-18170-66420-16180-64430-14190-62440-12200-60450-10210-58460-8220-56470-6230-54480-4240-52490-25000小船航行路线（数值解与解析解）②10.5v修改M文件boat.mfunctiondx=boat(t,x)v1=0.5;v2=2;dx=[v1-v2*x(1)/sqrt(x(1)^2+x(2)^2);-v2*x(2)/sqrt(x(1)^2+x(2)^2)];编写程序ts=0:60;x0=[0,-100];[t,x]=ode45(@boat,ts,x0);yt=-100:0;xt=50*((yt/(-100)).^0.75-(yt/(-100)).^1.25);plot(xt,yt,'b',x(1:55,1),x(1:55,2),'r')当54ts时，68.55100x，63.64100y，可以认为此时小船已到达B点，则渡河时间54ts。任意时刻小船位置txytxy00.0000-100.0000278.9700-46.271110.4950-98.0000289.0790-44.309720.9797-96.0001299.1668-42.352631.4541-94.0002309.2326-40.400341.9178-92.0006319.2754-38.453352.3706-90.0011329.2943-36.512062.8124-88.0020339.2883-34.577173.2428-86.0032349.2562-32.649283.6616-84.0048359.1969-30.729094.0685-82.0070369.1089-28.8174104.4633-80.0098378.9907-26.9154114.8456-78.0133388.8409-25.0239125.2152