How to Repeat an Iteration in a For loop, Python 3

Question

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.

It's just `i = 1` / `while i <= n:` then, at the end, do `i += 1` ONLY IF you want to bump the value. — Tim Roberts, Apr 04 '22 at 23:54
@TimRoberts - I am not sure I understand sir, do you mind further illustrating your point using what I have so far? Thank you for your time. — , Apr 05 '22 at 00:01

score 0 · Accepted Answer · answered Apr 05 '22 at 00:30

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

So this looks right, however, I tried your suggestion and I got a 'ZeroDivsionError' for the line with variable 'er'. Any idea what could be wrong? — , Apr 05 '22 at 00:37
You may print(h) just before the 'er' computation. May be h become lower and lower and finish to zero ? — Olivier Pellier-Cuit, Apr 05 '22 at 00:41
I see that after applying print(h)... now i'm not sure how to mitigate this. — , Apr 05 '22 at 00:43
If h become lower and lower then chance that er will become greater and greater (because of 1/h) term. It's certainly not the excepted result. Perhaps check that K1, K2 is correctly computed (with your reference algorithm). May be there is an error somewhere — Olivier Pellier-Cuit, Apr 05 '22 at 00:58

Lutz Lehmann · Answer 2 · 2022-12-15T15:51:27.303

The idea using a dynamic loop is correct, you do not know how many steps are needed to pass the end of the integration interval.

With the same reason it makes no sense to pass or define an n at all.

You can twist the logic of the step size adaptation. You do not compute an optimal h, you just decide if the current h is acceptable. The step size update happens anyway at the end of each loop in preparation of the next method step computation. But only if the local unit-step error is smaller than the tolerance do you append the computed y value to the output list and advance the time.

A skeleton with some parts missing could look like this

def RKF45(f,t,y,t_final,TOL):
    T = [t]
    Y = [y]
    h = TOL**0.2
    while t < t_final:
        v4, v5, err = RKF45Step(f,t,y,h)
        err = 1e-20+abs(err)
        # scale the error by the expected solution magnitude
        if 0.01*TOL < 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
        q = max(0.1,min(2.5,q))
        h = q*h
        # min/max on h
    return np.asarray(T), np.asarray(Y)

How to Repeat an Iteration in a For loop, Python 3

2 Answers2