I am currently working on implementing the numerical method RKF45 (Runge-Kutta-Fehlberg-45) with adaptive step size into python 3 and I believe I am running into a fundamental loop issue that I cannot resolve. Note, the portion of this numerical method I am having trouble implementing is the adaptive step size. I understand the basic algorithm for how this may be implemented, which I will provide, but first lets take a look at the function I have constructed that performs the RF45 calculations:
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
Note that gf.f is a function I defined on a separate module that is given by:
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.
Advertisement
Answer
If i understand what you need:
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