当 scipy.integrate.odeint 中的输出达到 0 时停止积分答案

【问题标题】：Stop Integrating when Output Reaches 0 in scipy.integrate.odeint当 scipy.integrate.odeint 中的输出达到 0 时停止积分
【发布时间】：2019-05-11 08:31:59
【问题描述】：

我已经编写了一个代码，它可以查看物体在拖动时的抛射运动。我正在使用 scipy 中的 odeint 来执行前向欧拉方法。积分将一直运行，直到达到时间限制。当达到此限制或 ry = 0 的输出（即弹丸已着陆）时，我想停止积分。

def f(y, t):
    # takes vector y(t) and time t and returns function y_dot = F(t,y) as a vector
    # y(t) = [vx(t), vy(t), rx(t), ry(t)]

    # magnitude of velocity calculated
    v = np.sqrt(y[0] ** 2 + y[1] **2)

    # define new drag constant k
    k = cd * rho * v * A / (2 * m)

    return [-k * y[0], -k * y[1] - g, y[0], y[1]] 


def solve_f(v0, ang, h0, tstep, tmax=1000):
    # uses the forward Euler time integration method to solve the ODE using f function
    # create vector for time steps based on args
    t = np.linspace(0, tmax, tmax / tstep)

    # convert angle from degrees to radians
    ang = ang * np.pi / 180

    return odeint(f, [v0 * np.cos(ang), v0 * np.sin(ang), 0, h0], t)

solution = solve_f(v0, ang, h0, tstep)

我已经尝试了几个循环和类似的尝试在 ry = 0 时停止积分。并在下面找到了这个问题，但无法用 odeint 实现类似的东西。 solution[:,3] 是 ry 的输出列。有没有一种简单的方法可以用 odeint 做到这一点？

Does scipy.integrate.ode.set_solout work?

【问题讨论】：

标签： python scipy ode odeint

【解决方案1】：

结帐scipy.integrate.odehere。它比odeint 更灵活，可以帮助您完成您想做的事情。

一个使用垂直拍摄的简单示例，在接触地面之前一直集成：

from scipy.integrate import ode, odeint
import scipy.constants as SPC

def f(t, y):
    return [y[1], -SPC.g]

v0 = 10
y0 = 0

r = ode(f)
r.set_initial_value([y0, v0], 0)

dt = 0.1
while r.successful() and r.y[0] >= 0:
    print('time: {:.3f}, y: {:.3f}, vy: {:.3f}'.format(r.t + dt, *r.integrate(r.t + dt)))

每次调用r.integrate，r 都会存储当前时间和 y 值。如果要存储它们，可以将它们传递给列表。

【讨论】：

你可以对 odeint 做同样的事情，y = odeint(f, y, [t, t+dt])[-1] 执行大小为 dt 的时间步。
可以，但我发现面向对象的 ODE 集成器更适合控制集成。在您的解决方案中，您还需要使用最后一个解决方案（ode 自动执行）在每个步骤中指定初始条件
感谢使用 ode 的建议。我已经重写了我的函数，现在我可以在输出达到某个值时停止积分。但是集成似乎没有正确评估，我已经用修改后的代码更新了这个问题。谢谢，杰克。
别担心，我的角度是度数，现在全部排序！

【解决方案2】：

让我们将其作为边界值问题来解决。我们有条件x(0)=0, y(0)=h0, vx(0)=0, vy(0)=0 和y(T)=0。要有固定的积分区间，设置t=T*s，表示dx/ds=T*dx/dt=T*vx等。

def fall_ode(t,u,p):
    vx, vy, rx, ry = u
    T = p[0]
    # magnitude of velocity calculated
    v = np.hypot(vx, vy)
    # define new drag constant k
    k = cd * rho * v * A / (2 * m)

    return np.array([-k * vx, -k * vy - g, vx, vy])*T 

def solve_fall(v0, ang, h0):
    # convert angle from degrees to radians
    ang = ang * np.pi / 180
    vx0, vy0 = v0*np.cos(ang), v0*np.sin(ang)

    def fall_bc(y0, y1, p): return [ y0[0]-vx0, y0[1]-vy0, y0[2], y0[3]-h0, y1[3] ]

    t_init = [0,1]
    u_init = [[0,0],[0,0],[0,0], [h0,0]]
    p_init = [1]
    res = solve_bvp(fall_ode, fall_bc, t_init, u_init, p_init)
    print res.message
    if res.success: 
        print "time to ground: ", res.p[0]
    # res.sol is a dense output, evaluation interpolates the computed values
    return res.sol

sol = solve_fall(300, 30, 20)
s = np.linspace(0,1,501)
u = sol(s)
vx, vy, rx, ry = u
plt.plot(rx, ry)

【讨论】：