Прекратите интеграцию, когда выход достигнет 0 в scipy.integrate.odeint - PullRequest
Новая
       108

Прекратите интеграцию, когда выход достигнет 0 в scipy.integrate.odeint

0 голосов
/ 10 декабря 2018

Я написал код, который смотрит на движение снаряда с помощью перетаскивания.Я использую odeint от scipy, чтобы сделать прямой метод Эйлера.Интеграция выполняется до тех пор, пока не будет достигнут срок.Я хотел бы прекратить интеграцию, либо когда будет достигнут этот предел, или когда выход для 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?

Ответы [ 2 ]

0 голосов
/ 11 декабря 2018

Давайте решим это как проблему граничного значения.У нас есть условия 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)
0 голосов
/ 10 декабря 2018

Оформить заказ scipy.integrate.ode здесь .Он более гибкий, чем 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.Вы можете передать их в список, если хотите сохранить их.

...