Блог

Я изо всех сил пытаюсь понять, почему этот цикл while (ниже в коде) завершается только после 2 итераций.

Это из метода RungeKutta45 Fehlberg с размером шага adative (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.6543612251060553 e-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 цикл-буквально останавливается после 2-й итерации.

Вы понимаете, почему это так? Время 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
 

Спасибо!

На второй итерации попробуйте напечатать dt_new перед строкой t = 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 на третьей итерации больше не будет выполняться, что приведет к завершению.