python - 为什么这个 Python while 循环仅在 2 次迭代后终止?
问题描述
我正在努力理解为什么这个 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
谢谢!
解决方案
在第二次迭代中,尝试在此块中的 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 条件将不再成立,从而导致终止。
推荐阅读
- r - 如何缩小汇总预测数据和每月真实数据之间的差距。
- swift - Swift4:Alamofire 使用预签名链接上传到 s3 并检查 md5
- asp.net-mvc - 在 asp.net MVC 中使用 ValidateAntiForgeryToken 和 JQGrid?
- python - Django时间重新格式化
- python - 如何在没有函数的情况下将 executor.map 应用于 for 循环?
- google-apps-script - 使用谷歌应用脚本将生成的谷歌表单移动到特定的谷歌驱动器文件夹
- javascript - 在前端和后端使用常量整数
- sql-server - 如果有人对数据库执行任何查询,则发送电子邮件
- python - aiortc 集成到 Django 项目中
- node.js - 如何使用云功能在 Firestore 中输入数据?