为什么这个 Python while 循环仅在 2 次迭代后终止？

Question

我很难理解为什么这个 while 循环（在代码下方）仅在 2 次迭代后终止。

它来自自适应步长 RungeKutta45 Fehlberg 方法 (https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ss2017/numerische_Methoden_fuer_komplexe_Systeme_II/rkm-1.pdf) 第 10/11 页。

下面的结果是这个输出：

$ python3 runge_kutta_45_adaptive_optimalstepsizes.py 

这一步的错误：0.0

此步骤出错：1.6543612251060553e-24

while 循环的迭代次数为：2

上次 t 为：0.001

import numpy as np
import os
import matplotlib
from matplotlib import pyplot as plt

# using current y_n and t_n, finds the largest possible dt_new such that the TRUNCATION ERROR
# after 1 integration step with timestep = this new dt_new (the one we want to settle upon)
# REMAINS below some given desired accuracy epsilon_0 
# ADAPTIVE STEP-SIZE CONTROL
# always change dt to dt = h_new (optimal timestep change)

rhs_of_diff_Eq_str = "3 * t ** 2"

def first_derivative(t, y): # the first derivative of the function y(t)
    first_derivative_value = 3 * t ** 2
    return first_derivative_value

def get_RKF4_approx(t, y, dt):
    k1 = first_derivative(t,                  y                                                                           )
    k2 = first_derivative(t +  dt/4. ,        y +  dt*( (1./4.)*k1                                                 )      )
    k3 = first_derivative(t +  dt*(3./8.) ,   y +  dt*( (3./32.)*k1 + (9./32.)*k2                                  )      )
    k4 = first_derivative(t +  dt*(12./13.) , y +  dt*( (1932./2197.)*k1 - (7200./2197.)*k2 + (7296./2197.)*k3     )      )
    k5 = first_derivative(t + dt,             y +  dt*( (439./216.)*k1 - 8.*k2 + (3680./513.)*k3 - (845./4104)*k4  )      )
    RKF4 = y + dt * ( (25./216)*k1 + (1408/2565.)*k3 + (2197./4104.)*k4 - (1./5.)*k5  )
    return np.array([RKF4, k1, k2, k3, k4, k5])

def get_RKF5_approx_efficiently(t, y, dt, ks): # efficient ! re-uses derivative evaluations from RKF4 (previous) calculation.
    # ks is a numpy array
    # ks[0] is K1, ks[1] is K2, ... , ks[4] is K5
    k6 = first_derivative(t + dt*(1./2.),     y +  dt*(-(8./27.)*ks[0] + 2.*ks[1] - (3544./2565.)*ks[2] + (1859./4104.)*ks[3] - (11./40.)*ks[4]) )
    RKF5 = y + dt * (  (16./135.)*ks[0] + (6656./12825.)*ks[2] + (28561./56430.)*ks[3] - (9./50.)*ks[4]  +(2./55.)*k6  )
    return RKF5 # a number

ts = []
ys = []
tfinal = 10.0
nmax = 10**6
epsilon_0 = 10**(-6)
contor = 0
dt = 0.001
beta = 0.9
t = 0.0 # initial condition
y = 0.0 # initial condition

while (t < tfinal and contor < nmax): 
    contor += 1
    container_from_RKF4method = get_RKF4_approx(t, y, dt)
    RKF4 = container_from_RKF4method[0] # the RKF4 method's approximation for y_{n+1}
    ks =  container_from_RKF4method[1:]
    RKF5 = get_RKF5_approx_efficiently(t, y, dt, ks)
    error_at_this_step = abs(RKF5 - RKF4)
    print("error at this step: {}".format(error_at_this_step))

    if (error_at_this_step < epsilon_0 and error_at_this_step != 0.0):
        # yes, step accepted! need optimal timestep
        dt_new = beta * dt * (epsilon_0/error_at_this_step)**(0.25)
        ts.append(t)
        t += dt_new
        dt = dt_new
        y_new = RKF5
        ys.append(y_new)
        y = y_new
    else:
        if (error_at_this_step == 0.0): # it's perfect! keep carrying on with this timestep which gives 0 error.
            ts.append(t)
            t += dt
            y_new = RKF5
            ys.append(y_new)
            y = y_new
        else: # means that error_at_this_step > epsilon_0 and that error_at_this_step != 0
            # no, step not accepted. reiterate step using a lower timestep
            dt_new = beta * dt * (epsilon_0/error_at_this_step)**(0.2)
            dt = dt_new
            # no changes made to time t and y
            # repeat this step (reiterate step)
            # HERE THE PROBLEM SHALL BE! I DON'T KNOW WHY THE ABOVE 2 instructions are bad!

print("no of iterations of the while loop was: {}".format(contor))
ts = np.array(ts)
print("last time t was: {}".format(ts[-1]))

ys = np.array(ys)
plt.figure()
plt.plot(ts, ys, label='y values', color='red')
plt.xlabel('t')
plt.ylabel('y')
plt.title("RK45 adaptive step-size (optimal step-size always chosen) integration for dy/dt = f(y,t) \n" + "f(y,t) = " + rhs_of_diff_Eq_str)
plt.savefig("RK45_adaptive_step_size_optimal_step_size_results.pdf", bbox_inches='tight')

我尝试查看使用breakpoint() 并按下n 和/或s 的指令执行情况。 似乎while-loop 在第二次迭代后实际上停止了。

你明白为什么会这样吗？ 时间t 没有达到tfinal 或contor=10**6-1=nmax！

显示问题的pdb调试位是：

> /mnt/c/Users/iusti/Desktop/runge_kutta_45_adaptive_optimalstepsizes.py(46)<module>()
-> while (t < tfinal and contor < nmax):
(Pdb) s
> /mnt/c/Users/iusti/Desktop/runge_kutta_45_adaptive_optimalstepsizes.py(79)<module>()
-> print("no of iterations of the while loop was: {}".format(contor))
(Pdb) s
[2]+  Stopped                 python3 runge_kutta_45_adaptive_optimalstepsizes.py

谢谢！

Answer 1

在第二次迭代中，尝试在此块中的 t += dt_new 行之前打印 dt_new：

    if (error_at_this_step < epsilon_0 and error_at_this_step != 0.0):
    # yes, step accepted! need optimal timestep
    dt_new = beta * dt * (epsilon_0/error_at_this_step)**(0.25)
    ts.append(t)
    t += dt_new
    dt = dt_new
    y_new = RKF5
    ys.append(y_new)
    y = y_new

我猜 dt_new 的值会很大，将其添加到 t 会导致 t >= tfinal，因此第三次迭代的 while 条件将不再成立，从而导致终止。