以下のスクリプトは、閉じたパスの周りに磁力線を統合し、PythonでRunge-Kutta RK4を使用して元の値に戻ったときに停止します。 SciPy.integrate.odeint
を使用したいと思いますが、パスがほぼ閉鎖されているときにどのように停止するかわかりません。パスが閉じたときにSciPy.integrate.odeintを停止させるには?
もちろん、odeint
は、Pythonでの統合よりもはるかに高速かもしれませんが、私は盲目的に迂回して結果を閉じ込めることができますが、将来はもっと大きな問題を抱えていきます。
"OKが十分に近い - 今すぐ停止することができます。"メソッドをodeintに実装する方法はありますか?または、しばらくの間統合してチェックし、さらに統合し、チェックする必要があります。
discussionこのように、「SciPyからはできません」という答えがあるかもしれません。
注:私は通常、スピードアップのために所定のステッピングサイズでより正確なRK45(Runge-Kutta-Fehlberg)を使用しますが、ここでは簡単にしました。また、可変ステップサイズも可能にします。
更新:時々私は一定のステップサイズが必要です。私は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[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()
これは本当に素晴らしいですね!私はすぐにそれを試してみると、物事が行くことを知らせる。係数はよく見えますが、これは[RKF45](http://maths.cnam.fr/IMG/pdf/RungeKuttaFehlbergProof.pdf)ですか?私はいつもcythonを避けていましたが、これは導入するのに最適な方法です。素晴らしい答えをありがとう! – uhoh
うん、あなたはそれをrkf45/rkfと呼んでいます。 Cythonで幸運を祈る!明確でないものがある場合は教えてください。 – Julius
週末までにこれを試すことはできないと思うが、今は確かに受け入れることができる。あなたの助けをもう一度ありがとう! – uhoh