使用scipy.integrate.solve_bvp求解抛体运动遇网格节点超限问题

我是编程新手，尝试用scipy.integrate.solve_bvp求解受边界条件约束、考虑空气阻力的一维竖直抛体运动：已知发射角90°、落地时间，需要确定发射初速度和落地末速度。

我先通过忽略空气阻力的SUVAT公式得到初速度猜测值，编写了如下代码：

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

drag_coef = 0.47 # 高尔夫球平均阻力系数
air_density = 1.293
area = 1.45*10**-3 # 高尔夫球横截面积
mass = 45.9*10**-3 # 高尔夫球质量
g = -9.80665

def function(time, height):
    drag_factor = drag_coef * air_density * area / (2*mass)
    return height[1], (g - drag_factor*height[1]*np.abs(height[1]))

def boundary_conditions(height_0, height_end):
    return height_0[0], height_end[0]

time_scale = 10 # 落地时间
velocity_0_guess = 49 # 初速度猜测值

time = np.linspace(0, time_scale, time_scale*1000+1)
height_0 = np.zeros(len(time))
velocity_0 = velocity_0_guess * np.ones(len(time))
height = np.array((height_0, velocity_0))

res = integrate.solve_bvp(function, boundary_conditions, time, height, max_nodes = time_scale*1000+1)

print(res.y[1][0]) # 计算得到的初速度
print(res.y[1][time_scale*1000]) # 计算得到的落地末速度

print(res)

plt.plot(time, res.y[0], label="S_z")
plt.xlabel("time [s]")
plt.ylabel("displacement [m]")
plt.show()

plt.plot(time, res.y[1], label="V_z")
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
plt.show()

当time_scale=10、velocity_0_guess=49时，运行后返回success=False，错误信息为：

message: 'The maximum number of mesh nodes is exceeded.'

但将time_scale改为5、velocity_0_guess改为24.5时，算法可正常收敛，返回success=True。

我已查阅官方文档、Reddit、Stack Overflow并尝试使用ChatGPT，但问题仍未解决，现寻求技术帮助。

解决方案 #

1. 初始网格无需过密，让算法自适应调整

solve_bvp是自适应节点的边界值问题求解器，它会根据计算误差自动加密/稀疏网格。初始给过于密集的网格（比如10001个点）不仅没必要，还会让算法在迭代调整时容易触发max_nodes限制。建议初始用稀疏网格，比如10个点：

time = np.linspace(0, time_scale, 10)

2. 优化初始猜测值的分布

初始速度猜测不要用全常数，建议给一个更接近实际运动趋势的分布（比如从初速度猜测值线性过渡到落地时的负速度），帮助算法更快收敛：

# 初始速度从49线性降到-40（模拟落地时向下的速度）
velocity_0 = np.linspace(velocity_0_guess, -40, len(time))

3. 适当增大max_nodes参数

算法在自适应加密节点时可能需要更多节点，可适当调高max_nodes的值，比如设为20000：

res = integrate.solve_bvp(function, boundary_conditions, time, height, max_nodes=20000)

修改后的完整代码 #

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

drag_coef = 0.47
air_density = 1.293
area = 1.45*10**-3
mass = 45.9*10**-3
g = -9.80665

def function(time, height):
    drag_factor = drag_coef * air_density * area / (2*mass)
    return height[1], (g - drag_factor*height[1]*np.abs(height[1]))

def boundary_conditions(height_0, height_end):
    return height_0[0], height_end[0]

time_scale = 10
velocity_0_guess = 49

# 初始稀疏网格
time = np.linspace(0, time_scale, 10)
height_0 = np.zeros(len(time))
# 优化初始速度猜测分布
velocity_0 = np.linspace(velocity_0_guess, -40, len(time))
height = np.array((height_0, velocity_0))

# 增大max_nodes限制
res = integrate.solve_bvp(function, boundary_conditions, time, height, max_nodes=20000)

print("计算得到的初速度:", res.y[1][0])
print("计算得到的落地末速度:", res.y[1][-1])
print(res)

# 用求解结果的节点和值绘图
plt.plot(res.x, res.y[0], label="S_z")
plt.xlabel("time [s]")
plt.ylabel("displacement [m]")
plt.legend()
plt.show()

plt.plot(res.x, res.y[1], label="V_z")
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
plt.legend()
plt.show()

修改后运行，算法应该能正常收敛，返回success=True，同时得到正确的初速度和末速度。

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