Poniższy skrypt integruje linie pola magnetycznego wokół zamkniętych ścieżek i zatrzymuje się, gdy powraca do pierwotnej wartości w ramach pewnej tolerancji, używając Runge-Kutta RK4 w Pythonie. Chciałbym użyć SciPy.integrate.odeint, ale nie widzę, jak mogę go zatrzymać, gdy ścieżka jest w przybliżeniu zamknięta.Jak zmusić SciPy.integrate.odeint do zatrzymania po zamknięciu ścieżki?
Oczywiście, że odeint może być znacznie szybszy niż integracja w Pythonie, mógłbym pozwolić mu na luzie i szukać zamknięcia w wynikach, ale w przyszłości będę robić o wiele większe problemy.
Czy istnieje sposób, w jaki mogę zaimplementować "OK, który jest wystarczająco blisko - można zatrzymać teraz!" metoda w odint? Czy powinienem po prostu zintegrować na jakiś czas, sprawdzić, zintegrować więcej, sprawdzić ...
Ten discussion wydaje się istotne i wydaje się sugerować, że "nie można z wewnątrz SciPy" może być odpowiedź.
Uwaga: Zwykle używam RK45 (Runge-Kutta-Fehlberg), który jest bardziej dokładny przy danym rozmiarze stoku, aby go przyspieszyć, ale zachowałem to tutaj prosto. Umożliwia także zmianę wielkości kroku.
Aktualizacja: Ale czasami potrzebuję stałej wielkości kroku. Zauważyłem, że Scipy.integrate.ode zapewnia metodę testowania/zatrzymywania, ale wydaje się, że nie ma możliwości oceny w stałych punktach t. odeint umożliwia ocenę w stałych punktach t, ale wydaje się, że nie ma metody testowania/zatrzymywania.
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[i] )
k2 = F(x[i] + k1*h_over_two)
k3 = F(x[i] + k2*h_over_two)
k4 = F(x[i] + k3*h )
x[i+1] = x[i] + 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()
Hej, to wygląda naprawdę fantastycznie! Wkrótce spróbuję i powiem ci, że wszystko idzie. Współczynniki wyglądają znajomo, to jest [RKF45] (http://maths.cnam.fr/IMG/pdf/RungeKuttaFehlbergProof.pdf)? Zawsze udało mi się uniknąć cytonina, ale wygląda to na świetny sposób na przedstawienie. Dziękuję za wspaniałą odpowiedź !! – uhoh
Tak, to jest rkf45/rkf, jakkolwiek chcesz to nazwać. Powodzenia z Cythonem! Daj mi znać, jeśli jest coś, co nie jest jasne. – Julius
Nie sądzę, że będę mógł to wypróbować do weekendu, ale z pewnością mogę teraz zaakceptować. Jeszcze raz dziękuję za pomoc! – uhoh