Составляющей ода интеграции неизвестным пределом т, в зависимости от размера раствора

Я моделирую заряженные частицы, движущиеся через электромагнитное поле, и использую scipy ode. Код здесь упрощен, очевидно, но работает как пример. Проблема в том, что я хочу закончить интегрирование после ограничения на r, а не на t. поэтому интегрируйте dx/dt до точки, где норма (x) > r.

Я не хочу просто изменить функцию для интегрирования по r, однако, потому что позиция является функцией от t. могу ли я сделать определенный интеграл по несвязанной переменной или что-нибудь? 

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

def RHS(state, t, Efield, q, mp):

    ds0 = state[3]
    ds1 = state[4]
    ds2 = state[5]
    ds3 = q/mp * Efield*state[0]
    ds4 = q/mp * Efield*state[1]
    ds5 = q/mp * Efield*state[2]

    r = np.linalg.norm((state[0], state[1], state[2]))

    # if r > 30000 then do stop integration.....?

    # return the two state derivatives
    return [ds0, ds1, ds2, ds3, ds4, ds5]


ts = np.arange(0.0, 10.0, 0.1)
state0 = [1.0, 2.0, 3.0, 0.0, 0.0, 0.0]

Efield=1.0
q=1.0
mp=1.0

stateFinal = odeint(RHS, state0, ts, args=(Efield, q, mp))
print(np.linalg.norm(stateFinal[-1,0:2]))
1 2
python numerical-integration scipy

1 ответ:

Вы можете управлять процессом, шаг за шагом выполняя интеграцию с помощью scipy.интегрировать.ода . Пример: 

from scipy.integrate import ode
t_initial = 0
dt = 0.1
t_max = 10
r = ode(RHS).set_initial_value(state0, t_initial).set_f_params(Efield, q, mp)
solution = [state0]
while r.successful() and r.t < t_max:
    new_val = r.integrate(r.t + dt)
    solution.append(new_val)
    if np.linalg.norm((new_val[0], new_val[1], new_val[2])) > 30000:
        break
print(solution)
Обратите внимание, что сигнатура RHS должна быть изменена на def RHS(t, state, Efield, q, mp): для ode, независимая переменная идет первой, в отличие от odeint. Выход - это решение, вычисленное в приращениях dt независимой переменной, вплоть до момента окончания цикла (либо потому, что t_max достигнут, либо интегратор не работает, либо условие для break было выполнено). встречающийся).
1