数学建模笔记-斜抛运动建模

作者：Peiy_Liu
日期：2020-02-04

---

1. 数学模型

质量为m的质点作斜抛运动，假定质点出射后，只受到恒力$m\bf{g}$mg和空气阻力$\bf{f}$f作用，其中f的大小与速度大小的平方成正比，方向与速度相反，即$f=-kv^2 \tag{1.1}$f=−kv2(1.1)令$\bf{r}$r $= \left[\begin{matrix}x&y\end{matrix}\right]^T$=[x​y​]T为质点的位矢，由牛顿第二定律列出质点动力学方程${\bf{f}}+m{\bf{g}}=m\ddot{\bf r} \tag{1.2}$f+mg=mr¨(1.2)由于空气阻力方向与速度方向相反，$\bf{f}$f可写为${\bf{f}}=-kv^2{\bf{e}}=-k\left[\begin{matrix}\dot{x}\sqrt{\dot{x}^2+\dot{y}^2}\\\dot{y}\sqrt{\dot{x}^2+\dot{y}^2}\end{matrix}\right] \tag{1.3}$f=−kv2e=−k[x˙x˙2+y˙​2​y˙​x˙2+y˙​2​​](1.3)式中${\bf{e}}=\left[\begin{matrix}1/\sqrt{1+y'^2}&y'/\sqrt{1+y'^2}\end{matrix}\right]^T$e=[1/1+y′2​​y′/1+y′2​​]T为质点速度方向的单位矢量.由此，式（1.2）可写成一个二阶微分方程组$\begin{cases}m\ddot{x}=-k\dot{x}\sqrt{\dot{x}^2+\dot{y}^2}\\m\ddot{y}=-k\dot{y}\sqrt{\dot{x}^2+\dot{y}^2}-mg\end{cases}\tag{1.4}${mx¨=−kx˙x˙2+y˙​2​my¨​=−ky˙​x˙2+y˙​2​−mg​(1.4)再令$z_1=\dot{x},z_2=\dot{y}$z1​=x˙,z2​=y˙​方程组（1.4）化为形如${\bf Y'}={\bf f}(t,{\bf Y})\tag{1.5}$Y′=f(t,Y)(1.5)的标准形式$\begin{cases}\dot{z_1}=-kz_1\sqrt{z_1^2+z_2^2}/m\\\dot{z_2}=-kz_2\sqrt{z_1^2+z_2^2}/m-g\end{cases}\tag{1.6}${z1​˙​=−kz1​z12​+z22​​/mz2​˙​=−kz2​z12​+z22​​/m−g​(1.6)
就可以用MATLAB中提供的函数求解了。

2. MATLAB建模

2.1 MATLAB解常微分方程

很多常微分方程难以求出解析解，需要运用数值解法，如方程组（1.6）。MATLAB提供了如ode45、ode23、ode113等函数可以解出形如式（1.5）的常微分方程（组）的数值解，通用用法为

[T,Y] = solver(fun,tspan,y0)

solver可用ode45、ode23、ode113中的任意一个替代，fun为等式右边的$\bf f$f向量，它所代表的M文件须有如下形式

function dy = fun(t,y)
    dy = ...

不管是否在fun中用到，fun一定有两个参数t和y。tspan为求解区间向量，y0为初值向量。返回值中T是自变量的取值，Y的各列为数值解。对由一个高阶常微分方程化为的方程组。Y(:1)是对应t的初值解，Y(:n)为对应t处初值解的n-1阶导数。

2.2 编程求轨迹

先用上述方法解出$z_1,z_2$z1​,z2​（即$v_x,v_y$vx​,vy​），再对分速度分别进行数值积分则可得到轨迹上的点$(x,y)$(x,y),最后可使用plot()函数绘出轨迹图。下面是代码

m = 2;                                    %质点质量kg
k = 0.50;                                 %空气阻力f = -k * v^2，国际单位
g = 9.8;                                  %重力加速度
theta = pi / 4;                           %出射仰角rad
v0 = 20;                                  %初速度大小
vx0 = v0 * cos(theta);
vy0 = v0 * sin(theta);
z0 = [vx0, vy0];                          %构造微分方程组
x0 = 0;
y0 = 0;
f = @(t,z) [-k * z(1) * sqrt(z(1)^2 + z(2)^2) / m; -k * z(2) * sqrt(z(1)^2 + z(2)^2) / m - g];
[T, Z] = ode45(f, [0, 1.5], z0);
vx = Z(:,1);
vy = Z(:,2);
x = x0;
y = y0;
for i = 1:length(T)-1                     %手动数值积分
    tempx =  x(i) + vx(i) * (T(i+1) - T(i));
    tempy =  y(i) + vy(i) * (T(i+1) - T(i));
    x = [x;tempx];
    y = [y;tempy];
end
plot(x, y)

运行效果如图

##参考文献
司守奎《数学建模算法与应用》