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

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

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

---

1. 数学模型

质量为m的质点作斜抛运动，假定质点出射后，只受到恒力$m\mathbf{g}$和空气阻力$\mathbf{f}$作用，其中f的大小与速度大小的平方成正比，方向与速度相反，即$$f=-k{v}^2\text{\hspace{1em}(1.1)}$$则质点始终在初速度确定的竖直平面内运动，在质点运动平面内以高度为0处一点为原点，水平方向为$x$轴，竖直方向为$y$轴建立坐标系。令$\mathbf{r}$ $={\left[\begin{array}{cc}x& y\end{array}\right]}^T$为质点的位矢，由牛顿第二定律列出质点动力学方程$$\mathbf{f}+m\mathbf{g}=m\ddot{\mathbf{r}}\text{\hspace{1em}(1.2)}$$由于空气阻力方向与速度方向相反，$\mathbf{f}$可写为$$\mathbf{f}=-k{v}^2\mathbf{e}=-k\left[\begin{array}{c}\dot{x}\sqrt{{\dot{x}}^2+{\dot{y}}^2}\\ \dot{y}\sqrt{{\dot{x}}^2+{\dot{y}}^2}\end{array}\right]\text{\hspace{1em}(1.3)}$$式中$\mathbf{e}={\left[\begin{array}{cc}1/\sqrt{1+y^{\prime 2}}& y^{\prime }/\sqrt{1+y^{\prime 2}}\end{array}\right]}^T$为质点速度方向的单位矢量.由此，式（1.2）可写成一个二阶微分方程组$$\{\begin{array}{l}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{array}\text{\hspace{1em}(1.4)}$$再令$z_1=\dot{x},z_2=\dot{y}$方程组（1.4）化为形如$${\mathbf{Y}}^{\prime }=\mathbf{f}(t,\mathbf{Y})\text{\hspace{1em}(1.5)}$$的标准形式$$\{\begin{array}{l}\dot{z_1}=-k{z}_1\sqrt{z_1^2+z_2^2}/m\\ \dot{z_2}=-k{z}_2\sqrt{z_1^2+z_2^2}/m-g\end{array}\text{\hspace{1em}(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为等式右边的$\mathbf{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$（即$v_x,v_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)

运行效果如图

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