案例:蹦级运动员-符号运算求解

案例:蹦级运动员-符号运算求解

假设你受雇与某家蹦极公司。你的任务是要预测蹦极过程中在自由落体阶段蹦极运动员的速度，他说时间的函数。

得到的信息可以用更进一步的分析，如针对不同质量的蹦极运动员确定蹦极绳索的长度和必要强度。

基于牛顿第二定律和基本的物理、流体力学知识，建立数学模型:

其中  是集总阻力系数(kg/m)

%符号运算，对于常微分方程求解，dsolve函数
>> %注释，不执行
>> %help命令
>> help('dsolve')%英文状态下的单引号
%此处省略dsolve的用法...

>> syms m g cd v(t)
>> eq=diff(v,t)==g-cd/m*v^2;
>> sol=dsolve(eq)
 
sol =
 
                                            (g^(1/2)*m^(1/2))/cd^(1/2)//与t无关，舍去
                                           -(g^(1/2)*m^(1/2))/cd^(1/2)//与t无关，舍去
 -(g^(1/2)*m^(1/2)*tanh(cd^(1/2)*g^(1/2)*m^(1/2)*(C2 - t/m)))/cd^(1/2)

问题:一个质量为68.1kg的蹦极运动员从一个静止的热气球上滑落，计算蹦极运动员前12s自由落体过程中的速度。如果所使用的绳索无限长，那么还可以确定蹦极运动员能达到的最终速度。假设阻尼系数为0.25kg/m。

>> %初始条件：当t=0，v=0
>> cond=v(0)==0;
>> %方向键，向上的，可以快速输入以前输入的指令
>> sol=dsolve(eq,cond)
 
sol =
 
(g^(1/2)*m^(1/2)*tanh((cd^(1/2)*g^(1/2)*t)/m^(1/2)))/cd^(1/2)
 
>> m0=68.1;
>> g0=9.81;
>> cd0=0.25;
>> %符号表达式求值subs()
>> vt=subs(sol,'m',m0)
 
vt =
 
(10^(1/2)*681^(1/2)*g^(1/2)*tanh((10^(1/2)*681^(1/2)*cd^(1/2)*g^(1/2)*t)/681))/(10*cd^(1/2))
 

>> vt=subs(subs(subs(sol,'m',m0),'g',g0),'cd',cd0)
 
vt =
 
(10^(1/2)*681^(1/2)*981^(1/2)*tanh((10^(1/2)*681^(1/2)*981^(1/2)*t)/13620))/50

>> limS=double(limit(vt,'t',inf))

limS =

   51.6938

>> %命令行窗口中代码太多了，需要清屏，工作区的变量野太多了，也要清除
>> %clc，clear

>> %脚本文件的运行方法，直接输入脚本文件名称，不带后缀点

最后成型的M文件bengS.m  

%1.求解微分方程的通解

syms m g cd v(t)
eq=diff(v,t)==g-cd/m*v^2;%微分方程模型
cond=v(0)==0;%初值条件
sol=dsolve(eq,cond)%微分方程模型求解，通解

%2.符号表达式的求值

m0=68.1;
g0=9.81;
cd0=0.25;
vt=subs(subs(subs(sol,'m',m0),'g',g0),'cd',cd0) %把质量、g、阻力系数值带入
 i=1;
 
%3.循环求解在时间区间[0,20]内，自由落体速度值

 for ti=0:20
    val(i)=double(subs(vt,'t',ti));
    i=i+1;
 end

%4.对时间0-20s内的速度进行可视化绘图

tf=0:20;
%绘图函数plot
plot(tf,val,'-.','LineWidth',2)
title('自由落体速度随时间变化曲线')
xlabel('时间t(s)')
ylabel('速度(m/s)')
grid on
hold on
plot(12,val(13),'r.','MarkerSize',15)

%5,蹦极运动员的极限速度

limS=double(limit(vt,'t',inf))

 可视化图形

 补充说明

>> 对m文件里添加双重百分号%%可以对代码分块
>> 同时在界面左下角野会出现了脚本，对每一块做了说明