MATLAB2009_1_4降落伞

a2=g*(1-(v2/vt).^2);%第二个加速度
subplot(3,1,1)%选择第一个子窗口
plot(t,v1,'o',t,v2,'o')%画速度曲线
subplot(3,1,2)%选择第二个子窗口
plot(t,x1,'o',t,x2,'o')%画高度曲线
subplot(3,1,3)%选择第三个子窗口
(5)
其中
(6)
根据速度可求下落的高度。当kv< 1时，利用关系v= dx/dt，由(3)式可得
积分
即 。
当t= 0时，x= 0，可得常数为
降落伞下落的高度为
(7)
当初速度比较小时，就用双曲余弦之比的对数计算高度。
当kv> 1时，则由(5)式可得
积分
即
当t= 0时，x= 0，可得常数
降落伞下落的高度为
程序在执行时将调用一个函数zho源自8_2ode_fun1.m。%降落伞的加速度函数(2)
functionf=fun(t,v,flag,g,vt)
f=g-g*v^2/vt^2;%加速度
[说明](1)常微分方程数值解的ode45还可以带更多的输入参数，前面三个参数分别文件名、自变量和初始条件，第4个是标志，在这里没有用，第4个以后的参数是常数。这为求解初速度不同的微分方程的数值解带来了便利。
plot(t,a1,'o',t,a2,'o')%画加速度曲线
程序在执行时将调用一个函数zhou8_2ode_fun2.m。
%降落伞的速度和加速度函数
functionf=fun(t,x,flag,g,vt)
f=[x(2);%速度
g-g*x(2)^2/vt^2];%加速度
[说明]此程序显示了求二个一阶微分方程组的数值解的更一般的方法。显然，用微分方程的数值解计算这样的问题比用公式要简单很多。
gridon%加网格
ylabel('高度\itx\rm/m','fontsize',fs)%标记高度
subplot(3,1,3)%选择第三个子窗口
plot(t,a1,t,a2,'r')%画加速度曲线
gridon%加网格
xlabel('时间\itt\rm/s','fontsize',fs)%标记时间
ylabel('加速度\ita\rm/m\cdots^-^2','fontsize',fs)%标记加速度
alpha1=atanh(0);%第一个常数
alpha2=atanh(10/vt);%第二个常数
% v=inline('vt*coth(g*t/vt+a)','t','g','vt','a');%速度内线函数
% x=inline('vt^2/g*log(sinh(g*t/vt+a)/sinh(a))','t','g','vt','a');%高度内线函数
sv=diff(sx,'t');%求速度的符号解
sa=diff(sv,'t');%求加速度的符号解
x=subs(sx,{'g','v0','vt'},{g,0,vt});%将符号用常数代替
v=subs(sv,{'g','v0','vt'},{g,0,vt});%将符号用常数代替
a=subs(sa,{'g','v0','vt'},{g,0,vt});%将符号用常数代替
降落伞下降的规律
[问题]降落伞下降的规律
物体在空气中运动时，阻力的大小可以表示为F=CρAv2/2。其中ρ是空气的密度，A是物体的有效横截面积，C为阻力系数。一降落伞和人的极限速度为vT= 5m/s，当它从静止开始下落时，求它的速度和下落的高度。如果该降落伞开始没有打开，当速度达到v0= 10m/s时才打开，它的运动规律是什么？
a1=a(v1,g,vt);%第一个加速度
v2=v(t,g,vt,alpha2);%第二个速度
x2=x(t,g,vt,alpha2);%第二个高度
a2=a(v2,g,vt);%第二个加速度
figure%创建图形窗口
subplot(3,1,1)%选择第一个子窗口
plot(t,v1,t,v2,'r')%画速度曲线
[t2,X2]=ode45('zhou8_2ode_fun2',t,[0 10],[],g,vt);%求二阶微分方程的数值解
x1=X1(:,1);%第一个高度
v1=X1(:,2);%第一个速度
x2=X2(:,1);%第二个高度
v2=X2(:,2);%第二个速度
a1=g*(1-(v1/vt).^2);%第一个加速度
由(1*)式可得
牛顿第二定律是一个二阶微分方程，在这里化成两个一阶微分方程，以便于求数值解。
[程序]zhou8_2ode.m的第三部分如下。
%降落伞下落的高度,速度和加速度(求二阶微分方程的数值解)
[t1,X1]=ode45('zhou8_2ode_fun2',t,[0 0],[],g,vt);%求二阶微分方程的数值解
[程序]zhou8_2ode.m的第一部分如下。
%降落伞下落的高度,速度和加速度(用解析公式)
clear%清除变量
vt=5;%极限速度
g=9.8;%重力加速度
t=0:0.01:2;%时间向量
v=inline('vt*tanh(g*t/vt+a)','t','g','vt','a');%速度内线函数
x=inline('vt^2/g*log(cosh(g*t/vt+a)/cosh(a))','t','g','vt','a');%高度内线函数
(3)如P2_9之下图所示，初速度较小时，其加速度的方向与速度方向相同，并随着速度的增加而减少；初速度较大时，其加速度的方向与速度方向相反，大小也随着速度的增加而减少。不论初速度是大还是小，加速度最后都趋于0。
[算法]方法一：用解析式。由于MATLAB可计算复数，在画曲线时不计虚部，所以，不论初速度大小如何，速度都可以用(3)式或(5)式计算，高度既可以用(7)计算，也可以用(8)式计算。
[数学模型]
伞和人受到重力mg，方向向下；空气阻力F，方向向上。取向下为正方向，根据牛顿第二定律可列人和伞的运动方程
mg-Kv2=ma
其中K=CρA/2。由于a= dv/dt，可得微分方程
(1)
其中k2=K/mg。分离变量得
(2)
当kv< 1时，例如初速度为0的情况，积分上式
设当t= 0时，v=v0，可得常数为
(8)
当初速度比较大时，就用双曲正弦之比的对数计算高度。
[图示](1)如P2_9之上图所示，不论降落伞的初速度是小于极限速度还是大于极限速度，最后的速度都趋近于极限速度。
P2_9图
(2)如P2_9之中图所示，在初速度较大的情况下，在相同的时间内，降落伞下落的高度要大些，所以速度较大的高度曲线在速度较小的高度曲线上面。
利用反双曲正切函数 可得
速度为
由于 ，可知：当t→∞时，v→1/k，这就是极限速度vT。因此速度为
(3)
其中
(4)
当kv> 1时，例如初速度大于极限速度的情况，对(2)式积分
设当t= 0时，v=v0，可得常数为
利用反双曲正切函数可得
速度为
由于 ，可知：当t→∞时，v→1/k，这也是极限速度vT。因此速度为
v=subs(sv,{'g','v0','vt'},{g,10,vt});%将符号用常数代替
a=subs(sa,{'g','v0','vt'},{g,10,vt});%将符号用常数代替
x2=subs(x,'t',t);%高度
v2=subs(v,'t',t);%速度
a2=subs(a,'t',t);%加速度
alpha1=acoth(0);%第一个常数
alpha2=acoth(v20/vt);%第二个常数
两套指令计算的数据所画的曲线完全相同。
方法二：用一个一阶微分方程的数值解。由(1)式可得降落伞下落的加速度微分方程为
(1*)
这是一个关于速度的一阶微分方程，其初始条件为：当t= 0时v= 0或v= 10。该方程的数值解就是速度与时间的关系，通过积分可得高度与时间的关系。
[程序]zhou8_2ode.m的第二部分如下。
%降落伞下落的高度,速度和加速度(求一阶微分方程的数值解)
[t1,v1]=ode45('zhou8_2ode_fun1',t,0,[],g,vt);%求一阶微分方程的数值解(1)
[t2,v2]=ode45('zhou8_2ode_fun1',t,10,[],g,vt);%求一阶微分方程的数值解(1)
plot(t,v1,'.',t,v2,'.')%画速度曲线
subplot(3,1,2)%选择第二个子窗口
holdon%保持图像
plot(t,x1,'.',t,x2,'.')%画高度曲线