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

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.

CodePudding user response：

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