如何用Python实现炮弹运动建模，解决任务2的开发问题

已完成的任务1代码 #

任务1要求实现考虑空气阻力的炮弹运动微分方程，对应效果及代码如下：

import math
import numpy

# 补充原代码缺失的炮弹质量计算
r_cb = 0.15
rho_lead = 11343
mass_cb = rho_lead * (4/3) * math.pi * r_cb ** 3

def f(r, t):
    '''根据状态向量r和时间t，返回炮弹运动微分方程的导数'''
    x, y, vx, vy = r
    
    # 阻力系数、迎风面积、空气参数、重力加速度定义
    kappa = 0.47
    A = math.pi * r_cb ** 2
    rho_air = 1.23
    g = 9.81
    k = 0.5 * kappa * rho_air * A

    dx_dt = vx
    dy_dt = vy
    
    speed = numpy.hypot(vx, vy)
    dvx_dt = -k / mass_cb * speed * vx
    dvy_dt = -k / mass_cb * speed * vy - g
        
    return numpy.array([dx_dt, dy_dt, dvx_dt, dvy_dt])

任务2要求解析与实现 #

任务2要求计算不同发射角下的炮弹射程，找到最大射程对应的最优发射角，对应任务描述如下：

任务2核心要求 #

• 实现calculate_range(theta)函数：输入发射角theta（弧度），返回炮弹落地（y=0）时的水平射程
• 求解微分方程调用scipy.integrate.odeint，当计算到y首次小于0时停止，用线性插值计算精确落地位置
• 在20°到70°的发射角区间内扫描，找到最大射程对应的发射角，输出角度值（保留两位小数）和最大射程（保留两位小数）

任务2实现代码 #

from scipy.integrate import odeint
import numpy as np

def calculate_range(theta):
    v0 = 200  # 初始速度m/s
    # 初始状态向量：x=0, y=0, vx=v0*cos(theta), vy=v0*sin(theta)
    r0 = np.array([0, 0, v0 * np.cos(theta), v0 * np.sin(theta)])
    t_step = 0.01  # 时间步长，单位秒
    t = 0
    y_prev = 0
    x_prev = 0
    # 迭代计算直到y首次小于0
    while True:
        t += t_step
        sol = odeint(f, r0, [0, t])
        x_curr, y_curr = sol[-1][0], sol[-1][1]
        if y_curr < 0:
            # 线性插值计算精确落地时间和射程
            ratio = y_prev / (y_prev - y_curr)
            x_range = x_prev + ratio * (x_curr - x_prev)
            return x_range
        y_prev = y_curr
        x_prev = x_curr

# 扫描20°到70°的发射角，步长0.01度保证精度
theta_deg_list = np.arange(20, 70, 0.01)
max_range = 0
best_theta_deg = 0
for theta_deg in theta_deg_list:
    theta_rad = np.radians(theta_deg)
    current_range = calculate_range(theta_rad)
    if current_range > max_range:
        max_range = current_range
        best_theta_deg = theta_deg

# 输出结果
print(f"最优发射角：{best_theta_deg:.2f}°")
print(f"最大射程：{max_range:.2f}m")

内容的提问来源于stack exchange，提问作者Sean Ambrose