使用solve_ivp求解3D抛体运动ODE遇值错误，求解决方案

错误原因与修复 #

1. 落地判断索引错误 #

你在3D模型中仍沿用了原2D模型的落地判断索引：原2D模型中sol.y[2]对应z坐标，改为3D后z坐标的索引变为sol.y[4]（状态向量S为[x, vx, y, vy, z, vz]）。原代码中判断y坐标的符号变化，而y方向初始速度为0且无风时y坐标始终为0，导致np.where返回空数组，后续索引操作引发维度异常。

修复：将落地判断的索引从sol.y[2]改为sol.y[4]。

2. 函数调用参数缺失 #

same_distance_alfa函数需要传入B和c_windx两个必填参数，但你调用时仅传入B，导致参数传递异常。修复：调用时补充c_windx参数。

3. 距离计算逻辑修正 #

原代码仅返回x方向的落点距离，但3D情况下落点是(x,y,z)，z=0时需要同时记录x和y坐标，才能计算离散度。

修正后的完整代码 #

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 物理参数
m = 30
v0 = 301
g = 9.81
ro = 1.28
A = 0.0026
c = 1
k = 0.5 * c * ro * A
B = k / m
V = v0
c_max = 10

def dSdt(t, S, B, c_windx, c_windy):
    x, vx, y, vy, z, vz = S
    # 相对风速：物体速度与风速的差
    rel_vx = vx - c_windx
    rel_vy = vy - c_windy
    rel_vz = vz
    speed = np.hypot(np.hypot(rel_vx, rel_vy), rel_vz)
    return [
        vx,
        -B * speed * rel_vx,
        vy,
        -B * speed * rel_vy,
        vz,
        -B * speed * rel_vz - g
    ]

def get_landing_pos(angle, B, c_windx, c_windy, V=v0, t_max=600):
    """计算给定角度下的落点坐标(x,y)"""
    angle_rad = np.radians(angle)
    v0x = V * np.cos(angle_rad)
    v0y = 0  # 初始y方向速度为0
    v0z = V * np.sin(angle_rad)
    y0 = [0, v0x, 0, v0y, 0, v0z]
    
    sol = solve_ivp(
        dSdt, [0, t_max], y0, 
        t_eval=np.linspace(0, t_max, 10000),
        args=(B, c_windx, c_windy),
        atol=1e-7, rtol=1e-4
    )
    
    # 寻找z坐标从正变负的点（落地时刻）
    z_vals = sol.y[4]
    sign_changes = np.where(np.diff(np.sign(z_vals)) < 0)[0]
    if len(sign_changes) == 0:
        return None  # 未落地（时间不够）
    idx = sign_changes[0]
    # 线性插值计算精确落点
    t1, t2 = sol.t[idx], sol.t[idx+1]
    x1, x2 = sol.y[0][idx], sol.y[0][idx+1]
    y1, y2 = sol.y[2][idx], sol.y[2][idx+1]
    z1, z2 = z_vals[idx], z_vals[idx+1]
    # 计算z=0时的t
    t_land = t1 + (-z1)*(t2-t1)/(z2-z1)
    x_land = x1 + (t_land - t1)*(x2-x1)/(t2-t1)
    y_land = y1 + (t_land - t1)*(y2-y1)/(t2-t1)
    return (x_land, y_land)

def find_matching_angle(target_angle, B, c_windx, c_windy, tol=10, start_offset=5):
    """寻找与目标角度落点距离小于tol的另一个角度"""
    target_pos = get_landing_pos(target_angle, B, c_windx, c_windy)
    if target_pos is None:
        return None
    target_x, target_y = target_pos
    
    # 从target_angle+start_offset开始遍历角度
    for angle in np.arange(target_angle + start_offset, 90, 0.1):
        current_pos = get_landing_pos(angle, B, c_windx, c_windy)
        if current_pos is None:
            continue
        current_x, current_y = current_pos
        dist = np.hypot(current_x - target_x, current_y - target_y)
        if dist < tol:
            return angle
    return None

# ---------------------- 计算离散度 ----------------------
# 无风情况参数
windless_params = (B, 0, 0)
# 有风情况参数（示例：x方向风10m/s，y方向风5m/s）
windy_params = (B, 10, 5)

alfa = 15
# 无风时找匹配角度
beta_windless = find_matching_angle(alfa, *windless_params)
# 有风时找匹配角度
beta_windy = find_matching_angle(alfa, *windy_params)

if beta_windless is None or beta_windy is None:
    print("未找到匹配角度，请调整参数或时间范围")
else:
    # 获取无风时两个角度的落点
    pos_alfa_windless = get_landing_pos(alfa, *windless_params)
    pos_beta_windless = get_landing_pos(beta_windless, *windless_params)
    # 获取有风时两个角度的落点
    pos_alfa_windy = get_landing_pos(alfa, *windy_params)
    pos_beta_windy = get_landing_pos(beta_windy, *windy_params)
    
    # 计算离散度：落点之间的欧氏距离差（有风与无风的差值）
    disp_windless = np.hypot(
        pos_beta_windless[0] - pos_alfa_windless[0],
        pos_beta_windless[1] - pos_alfa_windless[1]
    )
    disp_windy = np.hypot(
        pos_beta_windy[0] - pos_alfa_windy[0],
        pos_beta_windy[1] - pos_alfa_windy[1]
    )
    disp_diff = disp_windy - disp_windless
    
    print(f"无风时，角度{alfa}°与{beta_windless:.1f}°的落点离散度：{disp_windless:.2f}m")
    print(f"有风时，角度{alfa}°与{beta_windy:.1f}°的落点离散度：{disp_windy:.2f}m")
    print(f"有风相对于无风的离散度变化：{disp_diff:.2f}m")

# ---------------------- 可视化角度-距离曲线（3D投影到x-y平面） ----------------------
angles = np.linspace(0, 90, 200)
# 无风时的x落点
x_windless = [get_landing_pos(ang, *windless_params)[0] for ang in angles if get_landing_pos(ang, *windless_params)]
# 有风时的x落点
x_windy = [get_landing_pos(ang, *windy_params)[0] for ang in angles if get_landing_pos(ang, *windy_params)]

plt.figure(figsize=(10,6))
plt.plot(angles[:len(x_windless)], x_windless, label='无风')
plt.plot(angles[:len(x_windy)], x_windy, label='有风')
plt.xlabel('发射角度(°)')
plt.ylabel('x方向落点距离(m)')
plt.title('发射角度与落点距离关系')
plt.legend()
plt.axvline(alfa, ls='--', color='r', label=f'目标角度{alfa}°')
plt.axvline(beta_windless, ls=':', color='r', label=f'无风匹配角度{beta_windless:.1f}°')
plt.axvline(beta_windy, ls=':', color='b', label=f'有风匹配角度{beta_windy:.1f}°')
plt.legend()
plt.show()

关键说明 #

1. 相对风速修正：原代码中相对风速的计算应为物体速度 - 风速（风对物体的阻力基于相对风速），修正了rel_vx = vx - c_windx等计算逻辑。
2. 落点精确计算：通过线性插值计算z=0时的精确x、y坐标，避免直接取相邻点平均的误差。
3. 离散度定义：这里的离散度指两个角度落点之间的欧氏距离，你可根据需求调整为x/y方向的偏差、或落点相对于无风时的偏移量等。
4. 参数灵活性：可修改windy_params中的风速值，测试不同风况下的离散度变化。

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