利用python实现数据可视化（二）：考虑空气阻力的抛体运动方程与轨迹方程可视化（Numpy库、Matplotlib库）-CSDN博客

一、考虑空气阻力

（一）二次方阻力模型： $f_r = -kv^2(k>0)$

初速度大小为 $v_0$ ，质量为 $m$ ，抛射角为 $\theta$ 的抛体，在二次方空气阻力的作用下，由牛顿第二定律可得：①在水平方向上，有 $ma_x = -kv_x^2$

②在垂直方向上，有 $ma_y = -mg-kv_y^2$

因为 $\frac{\mathrm{d} v}{\mathrm{d} t} = a$ ，

所以可得：③ $\frac{\mathrm{d} v_x}{v_x^2} =-k \frac{\mathrm{d} t}{m}$

④ $\frac{\mathrm{d} v_y}{mg+kv_y^2} = -\frac{\mathrm{d} t}{m}$

当时间t处于0时刻时，即t=0时，抛体在水平方向上的速度大小为初速度大小v0在水平方向上的分量，即 $v_0cos\theta$ ；在垂直方向上的速度大小为初速度大小v0在垂直方向上的分量，即 $v_0sin\theta$ ；

当时间t处于任意时刻t，即t=t时，抛体在水平方向上的速度大小为 $v_x$ ，在垂直方向上的速度大小为 $v_y$ ；

对③式取定积分，有 $\int_{v_0cos\theta }^{v_x} \frac{\mathrm{d} v_x}{ v_x^2} = (-\frac{k}{m})\int_{0}^{t} dt$

对④式取定积分，有 $\int_{v_0sin\theta }^{v_y}\frac{\mathrm{d}v_y }{ mg+kv_y^2} = -\frac{\int_{0}^{t}dt}{m}$

求解上述定积分，可得：⑤ $v_x = \frac{mv_0cos\theta }{m+ktv_0cos\theta }$

⑥ $v_y =\sqrt{ \frac{mg}{k}}\tan \left [ \arctan \left ( v_0sin\theta \sqrt{\frac{k}{mg}} \right )-t\sqrt{\frac{gk}{m}} \right ]$

又因为 $\frac{\mathrm{d} x}{\mathrm{d} t} = v_x$ 和 $\frac{\mathrm{d} y}{\mathrm{d} t} = v_y$ ，

所以有⑦： $dx= \frac{mv_0cos\theta }{m+ktv_0cos\theta }dt$

⑧： $dy=\sqrt{ \frac{mg}{k}}\tan \left [ \arctan \left ( v_0sin\theta \sqrt{\frac{k}{mg}} \right )-t\sqrt{\frac{gk}{m}} \right ]dt$

且当时间t等于0时刻时，即t=0时，抛体在水平方向上的位移大小为x=0，在垂直方向上的位移大小为y=0；当时间t等于任意时刻t时，即t=t时，抛体在水平方向上的位移大小为x=x，在垂直方向上的位移大小为y=y；

对⑦式取定积分，有 $\int_{0}^{x}dx=\int_{0}^{t} \frac{mv_0cos\theta }{m+ktv_0cos\theta }dt$

对⑧式取定积分，有

$\int_{0}^{y}dy=\int_{0}^{t}\sqrt{ \frac{mg}{k}}\tan \left [ \arctan \left ( v_0sin\theta \sqrt{\frac{k}{mg}} \right )-t\sqrt{\frac{gk}{m}} \right ]dt$

求解上述定积分，可得，

抛体的运动方程为： $x =\frac{m}{k}\ln \frac{m+ktv_0cos\theta }{m}$

$y=\frac{m}{k}\ln\frac{\cos \left [ arctan(v_0sin\theta\sqrt{\frac{k}{mg}} )-t\sqrt{\frac{gk}{m}} \right ]}{cos [arctan(v_0sin\theta \sqrt{\frac{k}{mg}})]}$

联立抛体的运动方程，消去时间t，即可得抛体的轨迹方程为：

$y=\frac{m}{k}\ln\frac{\cos \left [ arctan(v_0sin\theta\sqrt{\frac{k}{mg}} )-(e^\frac{kx}{m}-1)\sqrt{\frac{mg}{k}}/(v_0cos\theta) \right ]}{cos [arctan(v_0sin\theta \sqrt{\frac{k}{mg}})]}$

通过上述运算，我们可以得到二次阻力模型中的抛体的运动方程和轨迹方程，然后，通过编写相对应的python代码，就可以将抛体的运动方程和轨迹方程可视化，具体代码如下，可视化结果请见图1-1（运动方程可视化）和图1-2（轨迹方程可视化）

def quadratic_drag_model(m, g, k, v0, cita, t):
    """二次方阻力模型（f = -k * (v**2)）
       参数意义：
              m:物体的质量
              g:重力加速度
              k:阻力系数（k>0）
              v0:物体的初速度
              cita:物体的抛射角（弧度）
              t:物体运动的末时刻"""

    # 初始化存储列表
    list_t = numpy.linspace(0, t+1, 10000)  # 增加采样点，使曲线更平滑
    list_x = []  # 水平位移
    list_y = []  # 垂直位移
    list_vx = []  # 水平速度
    list_vy = []  # 垂直速度
    valid_t = []  # 有效时间点（物体落地前）

    # 预计算常数，减少重复计算
    cos_cita = numpy.cos(cita)
    sin_cita = numpy.sin(cita)
    sqrt_k_mg = numpy.sqrt(k / (m * g))  # 角度相关常数
    initial_angle = numpy.arctan(v0 * sin_cita * sqrt_k_mg)  # 初始角度
    sqrt_gk_m = numpy.sqrt(g * k / m)  # 角度随时间变化的系数
    sqrt_mg_k = numpy.sqrt(m * g / k)  # 垂直速度系数

    # 标记物体是否已落地
    has_landed = False

    for time in list_t:
        if has_landed:
            break  # 落地后停止计算

        # 1. 计算水平速度和位移（物理意义明确，无奇点）
        denominator = m + k * v0 * cos_cita * time
        if denominator <= 0:
            break  # 避免除以零或负值
        vx = m * v0 * cos_cita / denominator
        x = (m / k) * numpy.log(denominator / m)  # 水平位移

        # 2. 计算垂直方向角度（核心改进：限制角度范围，避免tan奇点）
        current_angle = initial_angle - time * sqrt_gk_m
        # 当角度接近π/2时（tan趋于无穷），认为物体已落地
        if current_angle >= numpy.pi / 2 - 1e-4:  # 保留微小余量避免数值爆炸
            has_landed = True
            break

        # 3. 计算垂直速度和位移
        vy = sqrt_mg_k * numpy.tan(current_angle)
        # 垂直位移（确保cos为正，避免log负数）
        cos_current = numpy.cos(current_angle)
        cos_initial = numpy.cos(initial_angle)
        if cos_current <= 0 or cos_initial <= 0:
            has_landed = True
            break
        y = (m / k) * numpy.log(cos_current / cos_initial)

        # 4. 物理约束：y必须非负（落地后停止）
        if y < -1e-6:  # 允许微小负值（数值误差）
            has_landed = True
            break

        # 保存有效数据
        valid_t.append(time)
        list_x.append(x)
        list_y.append(y)
        list_vx.append(vx)
        list_vy.append(vy)

    # 绘制四幅子图
    fig = matplotlib.pyplot.figure(figsize=(12,10))
    ax = fig.subplots(2,2)

    # 1. x-t图
    ax[0][0].plot(valid_t, list_x, color='red', label='x=x(t)')
    ax[0][0].set_title('x = x(t)')
    ax[0][0].set_xlabel('t')
    ax[0][0].set_ylabel('x')
    ax[0][0].grid()
    ax[0][0].legend()

    # 2. y-t图
    ax[0][1].plot(valid_t, list_y, color='blue', label='y=y(t)')
    ax[0][1].set_title('y = y(t)')
    ax[0][1].set_xlabel('t')
    ax[0][1].set_ylabel('y')
    ax[0][1].grid()
    ax[0][1].legend()

    # 3. Vx-t图
    ax[1][0].plot(valid_t, list_vx, color='orange', label=r'$V_x=V_x(t)$')
    ax[1][0].set_title(r'$V_x = V_x(t)$')
    ax[1][0].set_xlabel('t')
    ax[1][0].set_ylabel(r'$V_x$')
    ax[1][0].grid()
    ax[1][0].legend()

    # 4. Vy-t图
    ax[1][1].plot(valid_t, list_vy, color='green', label=r'$V_y=V_y(t)$')
    ax[1][1].set_title(r'$V_y = V_y(t)$')
    ax[1][1].set_xlabel('t')
    ax[1][1].set_ylabel(r'$V_y$')
    ax[1][1].grid()
    ax[1][1].legend()

    # 绘制轨迹图 y-x
    fig = matplotlib.pyplot.figure()
    ax = fig.subplots()
    ax.plot(list_x, list_y, color='purple', label='y=y(x)')
    ax.set_title(' y = y(x)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid()
    ax.legend()

    matplotlib.pyplot.show()

# 进行测试
quadratic_drag_model(1,9.8,2.5,10,numpy.pi/4,2)

图1-1

图1-2

（二）一次方阻力模型： $f_r = -kv(k>0)$

由牛顿第二定律可得，

抛体在水平方向上有①： $ma_x=-kv_x$

抛体在垂直方向上有②： $ma_y=-mg-kv_y$

因为 $\frac{\mathrm{d} v}{\mathrm{d} t} = a$ ，

所以可得③： $\frac{\mathrm{d} v_x}{v_x} = -\frac{kdt}{m}$

④： $\frac{\mathrm{d} v_y}{mg+kv_y}=-\frac{dt}{m}$

当时间t处于任意时刻t，即t=t时，抛体在水平方向上的速度大小为 $v_x$ ，在垂直方向上的速度大小为 $v_y$ ；

所以，对③式取定积分，有： $\int_{v_0cos\theta }^{v_x}\frac{\mathrm{d} v_x }{v_x} = -\frac{k\int_{0}^{t}dt}{m}$

对④式取定积分，有： $\int_{v_0sin\theta }^{v_y}\frac{\mathrm{d} v_y}{mg+kv_y} = -\frac{\int_{0}^{t}dt}{m}$

求解上述定积分，可得⑤： $v_x=\frac{v_0cos\theta }{e^\frac{kt}{m}}$

⑥： $v_y = \frac{\left ( \frac{mg}{k}+v_0sin\theta \right )}{e^\frac{kt}{m}} - \frac{mg}{k}$

又因为 $\frac{\mathrm{d} x}{\mathrm{d} t} = v_x$ 和 $\frac{\mathrm{d} y}{\mathrm{d} t} = v_y$ ，且当时间t等于0时刻时，即t=0时，抛体在水平方向上的位移大小为x=0，在垂直方向上的位移大小为y=0；当时间t等于任意时刻t时，即t=t时，抛体在水平方向上的位移大小为x=x，在垂直方向上的位移大小为y=y；

所以，对⑤⑥两式分别再取定积分并求解，可得：

抛体的运动方程为：

$x=-\frac{mv_0cos\theta \left (e^\frac{-kt}{m}-1 \right )}{k}$

$y=-\frac{m(\frac{mg}{k}+v_0sin\theta )(e^\frac{-kt}{m}-1)}{k}-\frac{mgt}{k}$

联立上述两式，消去t，可得抛体的轨迹方程为：

$y=x(\frac{mg}{kv_0cos\theta }+tan\theta )+\frac{m^2g}{k^2}\ln (1-\frac{kx}{mv_0cos\theta })$

具体代码如下，运动方程可视化结果请见图1-3，轨迹方程可视化结果请见图1-4

def primary_drag_model(m,g,k,v0,cita,t):
    """一次方阻力模型（f = -k*v ）
       参数意义：
              m:物体的质量
              g:物体所受的重力加速度
              k:阻力系数（k>0）
              v0:物体的初速度
              cita:物体的抛射角（弧度）
              t:物体运动的末时刻"""

    # 水平位移列表
    list_x = []
    # 垂直位移列表
    list_y = []
    # 水平方向速度列表
    list_vx = []
    # 垂直方向速度列表
    list_vy = []

    list_t = numpy.linspace(0,t+1,10000)

    # 计算每一时刻物体的水平位移量的大小、垂直位移量的大小、水平方向速度大小和垂直方向速度大小并将其添加至相对应的列表中
    for i in list_t:
        list_x.append(m*v0*numpy.cos(cita)*(1-1/numpy.exp(k*i/m))/k)
        list_y.append((m/k) * (1-1/numpy.exp(k*i/m)) * ((m*g + k*v0*numpy.sin(cita))/k) - m*g*i/k)
        list_vx.append(v0*numpy.cos(cita)*(1/numpy.exp(k*i/m)))
        list_vy.append((m*g/k + v0*numpy.sin(cita)) * (1/numpy.exp(k*i/m)) - m*g/k)

    fig = matplotlib.pyplot.figure(figsize=(12,10))
    ax = fig.subplots(2,2)

    # 绘制水平方向运动图像，即水平方向位移的大小关于时间的函数图像
    ax[0][0].plot(list_t, list_x, color='red',label='x=x(t)')
    ax[0][0].set_title('x=x(t)')
    ax[0][0].set_xlabel('t')
    ax[0][0].set_ylabel('x')
    ax[0][0].grid()
    ax[0][0].legend()

    # 绘制垂直方向运动图像，即垂直方向位移的大小关于时间的函数图像
    ax[0][1].plot(list_t, list_y, color='blue',label='y=y(t)')
    ax[0][1].set_title('y=y(t)')
    ax[0][1].set_xlabel('t')
    ax[0][1].set_ylabel("y")
    ax[0][1].grid()
    ax[0][1].legend()

    # 绘制水平方向速度的大小随时间变化的图像
    ax[1][0].plot(list_t, list_vx, color='orange',label=r'$V_x=V_x(t)$')
    ax[1][0].set_title(r'$V_x=V_x(t)$')
    ax[1][0].set_xlabel('t')
    ax[1][0].set_ylabel(r'$V_x$')
    ax[1][0].grid()
    ax[1][0].legend()

    # 绘制垂直方向速度的大小随时间变化的图像
    ax[1][1].plot(list_t, list_vy, color='green',label=r'$V_y=V_y(t)$')
    ax[1][1].set_title(r'$V_y=V_y(t)$')
    ax[1][1].set_xlabel('t')
    ax[1][1].set_ylabel(r'$V_y$')
    ax[1][1].grid()
    ax[1][1].legend()

    # 单独在一个图表中绘制物体的轨迹图像，即物体的垂直位移量的大小关于水平位移量的大小的函数图像
    list_Y = []
    refer_list_x = []
    for i in list_x:
        # 由于后续会涉及到对数运算，所以应当先判断真数是否大于0
        if 1-k*i/(m*v0*numpy.cos(cita)) > 0:
            Y = i * (numpy.tan(cita) + (m * g / (k * v0 * numpy.cos(cita)))) + numpy.log(
                1 - k * i / (m * v0 * numpy.cos(cita))) * m * m * g / (k * k)
            # 由于物体的轨迹是具有实际意义的，所以舍去负数部分
            if Y>=0:
                refer_list_x.append(i)
                list_Y.append(Y)
    fig = matplotlib.pyplot.figure()
    ax = fig.subplots()
    ax.plot(refer_list_x, list_Y, color='purple', label='y=y(x)')
    ax.set_title('y=y(x)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid()
    ax.legend()

    matplotlib.pyplot.show()

# 测试
primary_drag_model(1,9.8,5,35,numpy.pi/4,10)

图1-3

图1-4

二、不考虑空气阻力

此时，空气阻力为： $f_r=0$

同理，由牛顿第二定律可得，

抛体在水平方向上的受力情况为①： $ma_x = 0$

抛体在垂直方向上的受力情况为②： $ma_y = -mg$

因为 $\frac{\mathrm{d}v_x }{\mathrm{d} t} = a_x$ ； $\frac{\mathrm{d} v_y}{\mathrm{d} t} = a_y$

且当时间t处于0时刻时，即t=0时，抛体在水平方向上的速度大小为初速度大小v0在水平方向上的分量，即 $v_0cos\theta$ ；在垂直方向上的速度大小为初速度大小v0在垂直方向上的分量，即 $v_0sin\theta$ ；

当时间t处于任意时刻t，即t=t时，抛体在水平方向上的速度大小为 $v_x$ ，在垂直方向上的速度大小为 $v_y$ ；

所以，对①②两式取定积分并求解，

可得③： $v_x = v_0cos\theta$

④： $v_y = -gt+v_0sin\theta$

又因为 $\frac{\mathrm{d} x}{\mathrm{d} t}=v_x$ ； $\frac{\mathrm{d} y}{\mathrm{d} t}=v_y$ ，且当时间t等于0时刻时，即t=0时，抛体在水平方向上的位移大小为x=0，在垂直方向上的位移大小为y=0；当时间t等于任意时刻t时，即t=t时，抛体在水平方向上的位移大小为x=x，在垂直方向上的位移大小为y=y；

所以，对③④两式再取定积分并求解，可得

抛体的运动方程为： $x=tv_0cos\theta$

$y = tv_0sin\theta -\frac{gt^2}{2}$

联立上述两式，消去t，可得

抛体的轨迹方程为： $y=-\frac{gx^2}{2v_0^2cos^2\theta }+x\tan \theta$

具体代码如下，抛体的运动方程可视化请见图2-1，轨迹方程可视化请见图2-2

def resistance_free_model(m,g,k,v0,cita,t,):
    """无阻力模型
       参数意义：
              m:物体的质量
              g:物体所受的重力加速度
              k:阻力系数（k=0）
              v0:物体的初速度
              cita:物体的抛射角（弧度）
              t:物体运动的末时刻"""

    # 在无阻力模型中，阻力系数恒为0
    k = 0
    # 水平位移列表
    list_x = []
    # 垂直位移列表
    list_y = []
    # 水平方向速度列表
    list_vx = []
    # 垂直方向速度列表
    list_vy = []

    list_t = numpy.linspace(0,t+1,10000)

    # 计算每一时刻物体的水平位移量的大小、垂直位移量的大小、水平方向速度大小和垂直方向速度大小并将其添加至相对应的列表中
    for i in list_t:
        list_x.append(i*v0*numpy.cos(cita))
        list_y.append(i*v0*numpy.sin(cita) - g*i*i/2)
        list_vx.append(v0*numpy.cos(cita))
        list_vy.append(v0*numpy.sin(cita) - g*i)

    fig = matplotlib.pyplot.figure(figsize=(12,10))
    ax = fig.subplots(2,2)

    # 绘制水平方向运动图像，即水平方向位移的大小关于时间的函数图像
    ax[0][0].plot(list_t, list_x, color='red',label='x=x(t)')
    ax[0][0].set_title('x=x(t)')
    ax[0][0].set_xlabel('t')
    ax[0][0].set_ylabel('x')
    ax[0][0].grid()
    ax[0][0].legend()

    # 绘制垂直方向运动图像，即垂直方向位移的大小关于时间的函数图像
    ax[0][1].plot(list_t, list_y, color='blue',label='y=y(t)')
    ax[0][1].set_title('y=y(t)')
    ax[0][1].set_xlabel('t')
    ax[0][1].set_ylabel("y")
    ax[0][1].grid()
    ax[0][1].legend()

    # 绘制水平方向速度的大小随时间变化的图像
    ax[1][0].plot(list_t, list_vx, color='orange',label=r'$V_x=V_x(t)$')
    ax[1][0].set_title(r'$V_x=V_x(t)$')
    ax[1][0].set_xlabel('t')
    ax[1][0].set_ylabel(r'$V_x$')
    ax[1][0].grid()
    ax[1][0].legend()

    # 绘制垂直方向速度的大小随时间变化的图像
    ax[1][1].plot(list_t, list_vy, color='green',label=r'$V_y=V_y(t)$')
    ax[1][1].set_title(r'$V_y=V_y(t)$')
    ax[1][1].set_xlabel('t')
    ax[1][1].set_ylabel(r'$V_y$')
    ax[1][1].grid()
    ax[1][1].legend()

    # 单独在一个图表中绘制物体的轨迹图像，即物体的垂直位移量的大小关于水平位移量的大小的函数图像
    list_Y = []
    refer_list_x=[]
    # 从list_x中依次获取x的数值，并计算其相对应的Y的数值
    for i in list_x:
        Y = i * numpy.tan(cita) - g * i * i / (2 * v0 * v0 * numpy.cos(cita) * numpy.cos(cita))
        # 由于物体的轨迹是具有实际意义的，所以舍去负数部分
        if Y>=0:
            refer_list_x.append(i)
            list_Y.append(Y)

    fig = matplotlib.pyplot.figure()
    ax = fig.subplots()
    ax.plot(refer_list_x,list_Y,color='purple',label='y=y(x)')
    ax.set_title('y=y(x)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid()
    ax.legend()

    matplotlib.pyplot.show()

# 测试
resistance_free_model(1,9.8,10,35,numpy.pi/4,10)

图2-1

图2-2

三、上述三种情况的轨迹比较

在上面，我们讨论了抛体在三种不同的阻力的情况下的运动方程与轨迹方程求解思想和步骤，同样，为了更直观地观察阻力对抛体的影响，我们在这里使用代码将三种不同阻力情况下的轨迹图像绘制在一个图层中，具体代码如下，可视化结果见图3-1

def trajectory_comparison_chart(m,g,k,v0,cita,t):
    """无阻力、一次方阻力、二次方阻力模型轨迹对比图
           参数意义：
                  m:物体的质量
                  g:重力加速度
                  k:阻力系数（k>0）
                  v0:物体的初速度
                  cita:物体的抛射角（弧度）
                  t:物体运动的末时刻"""

    list_t = numpy.linspace(0, t + 1, 10000)
    # 用来存放无阻力模型的水平位移与垂直位移大小
    list_x1 = []
    refer_list_x1=[]
    list_y1 = []

    # 用来存放一次方阻力模型的水平位移与垂直位移大小
    list_x2 = []
    refer_list_x2=[]
    list_y2 = []

    # 用来存放二次方阻力模型的水平位移与垂直位移大小
    list_x3 = []
    refer_list_x3 = []
    list_y3 = []

    for i in list_t:
        # 根据时间分别计算出无阻力、一次方阻力、二次方阻力模型的水平位移量的大小
        list_x1.append(i * v0 * numpy.cos(cita))
        list_x2.append(m * v0 * numpy.cos(cita) * (1 - 1 / numpy.exp(k * i / m)) / k)
        if (m + k * v0 * numpy.cos(cita) * i) / m > 0:
            list_x3.append(m * numpy.log((m + k * v0 * numpy.cos(cita) * i) / m) / k)

    # 根据无阻力模型的水平位移量大小计算出垂直位移量大小
    for i in list_x1:
        y1 = i * numpy.tan(cita) - g * i * i / (2 * v0 * v0 * numpy.cos(cita) * numpy.cos(cita))
        # 由于物体的运动轨迹是具有实际意义的，所以舍去负数
        if y1>=0:
            refer_list_x1.append(i)
            list_y1.append(y1)

    # 根据一次方阻力模型的水平位移量大小计算出垂直位移量大小
    for i in list_x2:
        # 由于涉及到对数运算，所以应当先判断真数是否大于0
        if 1 - k * i / (m * v0 * numpy.cos(cita)) > 0:
            y2 = i * (numpy.tan(cita) + (m * g / (k * v0 * numpy.cos(cita)))) + numpy.log(
                1 - k * i / (m * v0 * numpy.cos(cita))) * m * m * g / (k * k)
            # 物体的运动轨迹是具有实际意义的，舍去负数
            if y2>=0:
                refer_list_x2.append(i)
                list_y2.append(y2)

    # 根据二次方阻力模型的水平位移量大小计算出垂直位移量大小
    has_landed = False   # 用来判断物体是否落地，False表示没有落地，True表示已经落地
    for i in list_x3:
        # 如果物体已经落地，跳出整个循环
        if has_landed == True:
            break
        if numpy.cos(
                numpy.arctan(v0 * numpy.sin(cita) * numpy.sqrt(k / (m * g))) - (numpy.exp(k * i / m) - 1) * numpy.sqrt(
                        g * m / k) / v0 * numpy.cos(cita)) / numpy.cos(
                numpy.arctan(v0 * numpy.sin(cita) * numpy.sqrt(k / (m * g)))) > 0:
            y3 = m * numpy.log(numpy.cos(
                numpy.arctan(v0 * numpy.sin(cita) * numpy.sqrt(k / (m * g))) - (numpy.exp(k * i / m) - 1) * numpy.sqrt(
                    g * m / k) / v0 * numpy.cos(cita)) / numpy.cos(
                numpy.arctan(v0 * numpy.sin(cita) * numpy.sqrt(k / (m * g))))) / k
            if y3 >= 0:
                refer_list_x3.append(i)
                list_y3.append(y3)
            else:
                # 如果垂直方向的位移量小于0，表示物体已经落地，后续的数值就不必计算了
                has_landed = True


    fig = matplotlib.pyplot.figure()
    ax = fig.subplots()
    ax.plot(refer_list_x1,list_y1,color='red',label=r'$f_r=0$')
    ax.plot(refer_list_x2,list_y2,color='blue',label=r'$f_r=-kv$')
    ax.plot(refer_list_x3,list_y3,color='purple',label=r'$f_r=-kv^2$')

    ax.set_title('Trajectory Comparison Chart')
    ax.set_xlabel('Horizontal Displacement')
    ax.set_ylabel('Vertical Displacement')
    ax.grid()
    ax.legend()

    matplotlib.pyplot.show()

# 测试
trajectory_comparison_chart(1,9.8,2.5,10,numpy.pi/4,2)

图3-1

根据图3-1，我们可以清晰地看出，在质量相同、所受重力加速度大小相同、初速度大小相同以及抛射角相同的情况下，不同空气阻力对抛体的运动轨迹的影响是很大的，在不考虑空气阻力的情况下，抛体的运动轨迹是一个理想的抛物线（图3-1中的红色曲线），而随着阻力的不断变化，抛体的运动轨迹也随之而变化，运动轨迹不再是一个理想的抛物线，而是逐渐变化为图3-1中的蓝色或紫色曲线，而后面的这些情况，也更加地与我们的实际生活相契合。

利用python实现数据可视化（二）：考虑空气阻力的抛体运动方程与轨迹方程可视化（Numpy库、Matplotlib库）

一、考虑空气阻力

（一）二次方阻力模型 ：

（二）一次方阻力模型：

二、不考虑空气阻力

三、上述三种情况的轨迹比较

（一）二次方阻力模型： $f_r = -kv^2(k>0)$

（二）一次方阻力模型： $f_r = -kv(k>0)$