我目前正在努力在 python 3 中實作具有自適應步長的數值方法 RKF45 (Runge-Kutta-Fehlberg-45)，我相信我遇到了一個我無法解決的基本回圈問題。請注意，我在實作此數值方法時遇到困難的部分是自適應步長。我了解如何實作它的基本演算法，我將提供，但首先讓我們看看我構建的執行 RF45 計算的函式：
def rkf45(n): # here we perform the necessary RKF45 computations
t0 = 0
t1 = 5 # this is just a label for the endpoint, not the i = 1 point.
y0 = 0
TOL = 5e-7
h = (t1 - t0) / n
vect = [0] * (n 1)
vectw = [0] * (n 1)
vect[0] = t = t0
vectw[0] = y = y0
for i in range(1, n 1):
k1 = h * gf.f(t, y)
k2 = h * gf.f(t (1/4) * h, y (1/4) * k1)
k3 = h * gf.f(t (3/8) * h, y (3/32) * k1 (9/32) * k2)
k4 = h * gf.f(t (12/13) * h, y (1932/2197) * k1 - (7200/2197) * k2 (7296/2197) * k3)
k5 = h * gf.f(t h, y (493/216) * k1 - 8 * k2 (3680/513) * k3 - (845/4104) * k4)
k6 = h * gf.f(t (1/2) * h, y - (8/27) * k1 2 * k2 - (3544/2565) * k3 (1859/4104) * k4 - (11/40) * k5)
er = (1/h) * ((1/360) * k1 - (128/4275) * k3 - (2197/7540) * k4 (1/50) * k5 (2/55) * k6)
# adaptive step size test goes here
vect[i] = t = t0 i * h
vectw[i] = y = y ((16/135) * k1 (6656/12825) * k3 (28561/56430) * k4 - (9/50) * k5 (2/55) * k6)
return vect, vectw
請注意，這gf.f是我在一個單獨的模塊上定義的一個函式，由：
def f(t, y):
a = -3 * t * y ** 2
b = 1 / (1 t ** 3)
return a b
Now, where I have commented # adaptive step size goes here is where my question comes in: I need to test whether abs(er) > TOL and if this is true, update the current step size h by h = h * q where q = (TOL / (2 * abs(er))) ** (1 / 4) and repeat the current iteration with this updated step size until abs(er) < TOL. From there, I need to use this updated h in the next iteration.
I have tried using a while loop to achieve this but I am definitely not implementing this correctly; probably because I am new and making a silly mistake. I have also tried using an if statement to test whether or not abs(er) > TOL and from there update h but I do not believe this enforces the for loop to repeat the current iteration with an updated h.
uj5u.com熱心網友回復：
如果我了解您的需求：
def rkf45(n): # here we perform the necessary RKF45 computations
t0 = 0
t1 = 5 # this is just a label for the endpoint, not the i = 1 point.
y0 = 0
TOL = 5e-7
h = (t1 - t0) / n
vect = [0] * (n 1)
vectw = [0] * (n 1)
vect[0] = t = t0
vectw[0] = y = y0
for i in range(1, n 1):
# Will loop until break ("see End of loop ?" comment)
while True:
k1 = h * gf.f(t, y)
k2 = h * gf.f(t (1/4) * h, y (1/4) * k1)
k3 = h * gf.f(t (3/8) * h, y (3/32) * k1 (9/32) * k2)
k4 = h * gf.f(t (12/13) * h, y (1932/2197) * k1 - (7200/2197) * k2 (7296/2197) * k3)
k5 = h * gf.f(t h, y (493/216) * k1 - 8 * k2 (3680/513) * k3 - (845/4104) * k4)
k6 = h * gf.f(t (1/2) * h, y - (8/27) * k1 2 * k2 - (3544/2565) * k3 (1859/4104) * k4 - (11/40) * k5)
er = (1/h) * ((1/360) * k1 - (128/4275) * k3 - (2197/7540) * k4 (1/50) * k5 (2/55) * k6)
# End of loop ?
if abs(er) <= TOL:
break
# update h before starting new iteration
h *= (TOL / (2 * abs(er))) ** (1 / 4)
# Continue
vect[i] = t = t0 i * h
vectw[i] = y = y ((16/135) * k1 (6656/12825) * k3 (28561/56430) * k4 - (9/50) * k5 (2/55) * k6)
return vect, vectw
uj5u.com熱心網友回復：
使用動態回圈的想法是正確的，您不知道需要多少步才能通過積磁區間的末尾。
出于同樣的原因，傳遞或定義一個根本沒有意義n。
您可以扭曲步長適應的邏輯。您不計算最佳值h，您只需決定當前h是否可以接受。步長更新無論如何都會在每個回圈結束時發生，以準備下一個方法步驟計算。但只有當區域單位步長誤差小于容差時，您才會將計算y值附加到輸出串列并提前時間。
缺少某些部分的骨架可能看起來像這樣
def RKF45(f,t,y,t_final,TOL):
T = [t]
Y = [y]
while t < t_final:
v4, err = RK45step(f,t,y,h)
# scale the error by the expected solution magnitude
if err<TOL:
t,y = t h, y h*v4
T.append(t)
Y.append(y)
q = 0.9*(TOL/err)**0.25
# min/max on q
h = q*h
# min/max on h
return np.asarray(T), np.asarray(Y)
轉載請註明出處，本文鏈接：https://www.uj5u.com/caozuo/456869.html
標籤：python for-loop while-loop runge-kutta
上一篇：向Set<String>添加回應