数学实验第二次作业――常微分方程数值求解

实验4常微分方程数值解实验目的：1.练习数值积分的计算；2.掌握用MATLAB软件求微分方程初值问题数值解的方法；3.通过实例学习用微分方程模型解决简化的实际问题；4.了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。实验内容：3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。解答如下：这是一个典型的牛顿第二定律问题，分析火箭受力情况；先规定向上受力为正数建立数学模型：A燃料未燃尽前，在任意时刻（t60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv^2,k=0.4,v表示此时火箭速度；此时火箭收到的合力为F1=（F-mg-f）；火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻火箭质量为m=1400-18*t根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))由此可利用龙格-库塔方法来实现，程序实现如下Function[dx]=rocket[t,x]%建立名为rocket的方程m=1400;k=0.4;r=-18;g=9.8;%给出题目提供的常数值dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];%以向量的形式建立方程[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t);%给出a的表达式End;ts=0:60;%根据题目给定燃烧率计算出燃料燃尽的时间，确定终点x0=[0,0];%输入x的初始值[t,x]=ode15s(@rocket,ts,x0);%调用ode15s计算[t,x];h=x(:,1);v=x(:,2);plot(t,x(:,1)),grid;%绘出火箭高度与时间的关系曲线title('h-t');xlabel('t/s');ylabel('h/m'),pause;plot(t,x(:,2)),grid;%绘出火箭速度与时间的曲线关系title('v-t');xlabel('t/s');ylabel('v/m/s'),pause;a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t);plot(t,a),grid;%绘出火箭加速度与时间的曲线关系title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause火箭高度随时间变化的曲线火箭速度随时间变化的曲线火箭加速度随时间变化的曲线数据过多，故截取部分如下第一列为时间，第二列为火箭高度，第三列为火箭速度519823.748258.96775210083.18259.92335310343.59260.86135410604.94261.7945510867.21262.72215611130.39263.64625711394.49264.56715811659.49265.4825911925.42266.38986012192.25267.2908由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2B燃料燃尽之后，与A类似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv^2,k=0.4,v表示此时火箭速度；此时火箭收到的合力为F2=（-mg-f）；火箭的初始质量为320kg，恒定根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320;程序实现如下function[dx]=rocket2(t,x)%建立以rocket2为名的函数dx=[x(2);-9.8-0.4.*x(2).^2/320];%以向量的形式建立方程ts=60:120;%给出初始时刻，估计终点时刻x0=[12190,267.26];%给出x初始值[t,x]=ode15s(@rocket2,ts,x0);%调用ode15s计算[t,x]plot(t,x(:,1)),grid;%绘出火箭高度随时间变化的曲线title('h-t');xlabel('t/s'),ylabel('h/m'),pause;plot(t,x(:,2)),grid;%绘出火箭速度随时间的变化曲线title('v-t');xlabel('t/s'),ylabel('v/m/s'),pause;v=x(:,2);a=-9.8-0.4*v.^2/320;%给出加速度的具体表达式plot(t,a),grid;%绘出火箭加速度随时间变化的曲线title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause得到的曲线图形如下火箭高度随时间的变化曲线从图中可以大致看出，最高点在13km左右，火箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格大致如下0.00641.28220.0093-0.00210.00651.29050.0074-0.00170.00661.29710.0059-0.00140.00671.30230.0046-0.00120.00681.30630.0034-0.00110.00691.30920.0023-0.00100.00701.31100.0013-0.00100.00711.31170.0003-0.00100.00721.3116-0.0007-0.00100.00731.3104-0.0017-0.0010从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s^2;本题总结：这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定，先用ode15s速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长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）,情况又如何建立数学模型：在任意时刻t，小船位于（x，y），此时速度为v，根据物理中路程与速度的关系，知路程的微分为速度v，由此中，小船在x,y方向上的速度分别为：X：dxdt=v1-v2*x√x2+y^2Y:dydt=-v2*y√x2+y^2初始条件为dxdt(t0)=v1*（1）dydt(t0)=v2*（2）现求其解析解,利用微积分知识*（1）÷*（2）得dxdy=v1∗√x2+y^2−v2y+xy;*(3)将*（3）右端带根号部分分子分母均÷y得dxdy=v1∗√(x/y)2+1−v2+xy;令x/y=p，得到dx/dy=dp/dy*y+p;dx/dy=p;分离变量有y^(-k)=c(√1+p2+p);代入初值可确定当t=0时，y=-d,x=0,p=x/y=0,C=(-d)^(-k)y^(-k)/d^(-k)=√1+x2y^2+x/y;根据隐函数与显函数的关系可得到解析解：x=-d2[(y−d)^(-k)-(yd)^(1+k)]已知d=100m,v1=1m/s,v2=2m/s;先利用龙格库塔方法求解渡河时间，及任意时刻小船的位置（x，y），及航行曲线，与解析解比较此时，k=0.5，d=100用MATLB编写源程序如下：functiondx=boat(t,x)%建立名为boat的函数d=-100;v1=1;v2=2;%给定常数值s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s];%用向量形式建立方程ts=0:70;%大致估算，确定终点值，给定步长为”1”x0=[0,-100];%给出x初始值[t,x]=ode23s(@boat,ts,x0);%调用ode23s计算[t,x];plot(t,x),grid,%绘出x(t),y(t)的函数曲线（图21）gtext('x(t)');gtext('y(t)');pause;plot(t,x(:,1)),%绘出x随时间变化的曲线（图22）grid;xlabel('t/s');ylabel('x/m');pause;plot(t,x(:,2)),grid,%绘出y随时间变化的曲线（图23）title('y-t');xlabel('t/s'),ylabel('y/m');图21图22图23得到数据如下00-10010.990037-9821.959741-96.000232.908828-94.000843.836822-92.00254.742902-90.004265.627337-88.007376.489358-86.011987.327888-84.0184606.537502-1.74536615.596041-1.26513624.637422-0.86066633.664544-0.53305642.680323-0.28315651.687817-0.11154660.690109-0.018567########0682.15E-080691.63E-080从表格中数据可知，在大约67s时y=0即船到达对岸目的地，为比较，先进行解析解的求解设计程序如下：functionx=f(y)k=0.5;x=-0.5.*(-0.01).^k.*y.^(k+1)+0.5.*(-0.01).^(-k).*y.^(-k+1);y=[0:-0.1:-100];fori=0:1:1000;x(:,i+1)=xy(-i/10);endplot(x,y);grid;gtext('x');gtext('y');由此可以看出，由数值解和解析解得到的x-y曲线相差不多，所以可以认为解析解正确改变水流速度v1，只要在原有程序基础上重新复制给v1=0,0.5,1.5,2,同时适当改变终点值即可现实现程序如下A:v1=0，d=-100;v1=0;v2=2;s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s]ts=0:60;x0=[0,-100];[t,x]=ode15s(@boat21,ts,x0);[t,x];plot(x(:,1),x(:,2)),grid,title('y-x');pause;plot(t,x(:,1)),grid;xlabel('t/s');ylabel('x/m');plot(t,x(:,1)),grid;title('x-t');xlabel('t/s'),ylabel('x/m'),pause;plot(t,x(:,2)),grid,title('y-t');xlabel('t/s'),ylabel('y/m');图形如下从此图形中我们可以看到，船并未偏离x=0的点，我们也可以从直观想象中的得到，当水速v1=0时，只要出发时，船头对准目标点，船将一直朝着直线向目的地行进400-20410-18420-16430-14440-12450-10460-8470-6480-4490-2500########从表格中数据我们也可以很清楚地看到路程与时间是成明显的线性关系的，这是与我们水速为0的必然结果，由此也可以验证我们模型基本正确现改变水速B:v1=0.5程序实现如下d=-100;v1=0.5;v2