Как получить scipy.integrate.odeint, чтобы остановиться, когда путь закрыт? ⇐ Python

[Anonymous]

Редактировать: < /strong> прошло пять лет, есть scipy.integrate.odeint < /code> научился останавливаться? Я хотел бы использовать scipy.integrate.odeint , но я не могу видеть, как я могу сказать, чтобы остановиться, когда путь приблизительно закрыт. /> Есть ли способ, которым я могу реализовать метод « ok, который достаточно близко - вы можете остановиться сейчас! < /em>» в Odeint? Или я должен просто интегрироваться на некоторое время, проверить, интегрировать больше, проверить ... < /p>

Это обсуждение кажется актуальным, и, кажется, предполагает, что «вы не можете из Scipy» может быть ответом. здесь. Это также делает возможным размер шага переменной. < /P>

Обновление: < /strong> Но иногда мне нужен фиксированный размер шага. Я обнаружил, что scipy.integrate.ode предоставляет метод тестирования/остановки ode.solout (t, y) , но, похоже, не может оценивать в фиксированных точках t . Odeint позволяет оценку в фиксированных точках t , но, похоже, не имеет метода тестирования/остановки.



def rk4Bds_stops(x, h, n, F, fclose=0.1):

h_over_two, h_over_six = h/2.0, h/6.0

watching = False
distance_max = 0.0
distance_old = -1.0

i = 0

while i < n and not (watching and greater):

k1 = F( x )
k2 = F( x + k1*h_over_two)
k3 = F( x + k2*h_over_two)
k4 = F( x + k3*h )

x[i+1] = x + h_over_six * (k1 + 2.*(k2 + k3) + k4)

distance = np.sqrt(((x[i+1] - x[0])**2).sum())
distance_max = max(distance, distance_max)
getting_closer = distance < distance_old

if getting_closer and distance < fclose*distance_max:
watching = True

greater = distance > distance_old
distance_old = distance

i += 1

return i

def get_BrBztanVec(rz):

Brz = np.zeros(2)

B_zero = 0.5 * i * mu0 / a
zz = rz[1] - h
alpha = rz[0] / a
beta = zz / a
gamma = zz / rz[0]

Q = ((1.0 + alpha)**2 + beta**2)
k = np.sqrt(4. * alpha / Q)

C1 = 1.0 / (pi * np.sqrt(Q))
C2 = gamma / (pi * np.sqrt(Q))
C3 = (1.0 - alpha**2 - beta**2) / (Q - 4.0*alpha)
C4 = (1.0 + alpha**2 + beta**2) / (Q - 4.0*alpha)

E, K = spe.ellipe(k**2), spe.ellipk(k**2)

Brz[0] += B_zero * C2 * (C4*E - K)
Brz[1] += B_zero * C1 * (C3*E + K)

Bmag = np.sqrt((Brz**2).sum())

return Brz/Bmag

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

pi = np.pi
mu0 = 4.0 * pi * 1.0E-07

i = 1.0 # amperes
a = 1.0 # meters
h = 0.0 # meters

ds = 0.04 # step distance (meters)

r_list, z_list, n_list = [], [], []
dr_list, dz_list = [], []

r_try = np.linspace(0.15, 0.95, 17)

x = np.zeros((1000, 2))

nsteps = 500

for rt in r_try:

x[:] = np.nan

x[0] = np.array([rt, 0.0])

n = rk4Bds_stops(x, ds, nsteps, get_BrBztanVec)

n_list.append(n)

r, z = x[:n+1].T.copy() # make a copy is necessary

dr, dz = r[1:] - r[:-1], z[1:] - z[:-1]
r_list.append(r)
z_list.append(z)
dr_list.append(dr)
dz_list.append(dz)

plt.figure(figsize=[14, 8])
fs = 20

plt.subplot(2,3,1)
for r in r_list:
plt.plot(r)
plt.title("r", fontsize=fs)

plt.subplot(2,3,2)
for z in z_list:
plt.plot(z)
plt.title("z", fontsize=fs)

plt.subplot(2,3,3)
for r, z in zip(r_list, z_list):
plt.plot(r, z)
plt.title("r, z", fontsize=fs)

plt.subplot(2,3,4)
for dr, dz in zip(dr_list, dz_list):
plt.plot(dr, dz)
plt.title("dr, dz", fontsize=fs)

plt.subplot(2, 3, 5)
plt.plot(n_list)
plt.title("n", fontsize=fs)

plt.show()


Подробнее здесь: https://stackoverflow.com/questions/330 ... -is-closed