如何实现自适应步长Runge-Kutta Cash-Karp?
How to implement adaptive step size Runge-Kutta Cash-Karp?
尝试实现自适应步长 Runge-Kutta Cash-Karp 但失败并出现此错误:
home/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in double_scalars from ipykernel import kernelapp as app
我试图解决的 ODE(并在下面的示例中使用,从高阶转换为一阶 ODE 系统)如下:
这是我正在使用的实现:
def rkck(f, x, y, h, tol):
#xn = x + h
err = 2 * tol
while (err > tol):
xn = x + h
k1 = h*f(x,y)
k2 = h*f(x+(1/5)*h,y+((1/5)*k1))
k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
yn4 = y + ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
yn5 = y + ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
err = yn4[-1]-yn5[-1]
if (err != 0):
h = 0.8 * h * (tol/err)**(1/float(5))
yn = yn4
return xn, yn
def integrate_sStepControl(f, t0, y0, tend, h, tol):
T = [t0]
Y = [y0]
t = t0
y = y0
while (t < tend):
h = min(h, tend-t)
t, y = rkck(f, t, y, h, tol)
T.append(t)
Y.append(y)
return np.array(T), np.array(Y)
def f_1(t,y):
return np.array([ y[1], -y[0]-(y[0])**3 ])
Y0_f1 = np.array([1.0,1.0])
# Execution
h = 0.05
tv, yv = integrate_sStepControl(f=f_1, t0=0.0, y0=Y0_f1, tend=100.0, h=h, tol=1.0E-05)
print("[ %20.15f, %20.15f]"%(yv[-1,0], yv[-1,1]) )
plt.plot(tv, yv)
得到上面的错误,但它被绘制出来了。不知道我在这里做错了什么:-/
编辑:
添加了对 err == 0
的检查
您需要实际传回计算出的新步长 h
以便它可以在第一步的主循环中使用。
误差应按标准计算。添加一些小数以避免被零除。
四阶方法的主要误差项是 C*h^5
。这必须与所需的本地错误 tol*h
进行比较。总的来说,这会导致计算最优 h
的第 4 个根。取 5 次方根可以提供某种抑制作用,但是对全局误差的影响并不是很直接。
def rkck(f, x, y, h, tol):
#xn = x + h
err = 2 * tol
while (err > tol):
xn = x + h
k1 = h*f(x,y)
k2 = h*f(x+(1/5)*h,y+((1/5)*k1))
k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
dy4 = ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
dy5 = ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
err = 1e-2*tol+max(abs(dy4-dy5))
# h = 0.95 * h * (tol/err)**(1/5)
h = 0.8 * h * (tol*h/err)**(1/4)
yn = y+dy4
return xn, yn, h
def integrate_sStepControl(f, t0, y0, tend, h, tol):
T = [t0]
Y = [y0]
t = t0
y = y0
while (t < tend):
h = min(h, tend-t)
t, y, h = rkck(f, t, y, h, tol)
T.append(t)
Y.append(y)
return np.array(T), np.array(Y)
因此,您的示例给出了解决方案、时间步长和 error/tol
的以下图
尝试实现自适应步长 Runge-Kutta Cash-Karp 但失败并出现此错误:
home/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in double_scalars from ipykernel import kernelapp as app
我试图解决的 ODE(并在下面的示例中使用,从高阶转换为一阶 ODE 系统)如下:
这是我正在使用的实现:
def rkck(f, x, y, h, tol):
#xn = x + h
err = 2 * tol
while (err > tol):
xn = x + h
k1 = h*f(x,y)
k2 = h*f(x+(1/5)*h,y+((1/5)*k1))
k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
yn4 = y + ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
yn5 = y + ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
err = yn4[-1]-yn5[-1]
if (err != 0):
h = 0.8 * h * (tol/err)**(1/float(5))
yn = yn4
return xn, yn
def integrate_sStepControl(f, t0, y0, tend, h, tol):
T = [t0]
Y = [y0]
t = t0
y = y0
while (t < tend):
h = min(h, tend-t)
t, y = rkck(f, t, y, h, tol)
T.append(t)
Y.append(y)
return np.array(T), np.array(Y)
def f_1(t,y):
return np.array([ y[1], -y[0]-(y[0])**3 ])
Y0_f1 = np.array([1.0,1.0])
# Execution
h = 0.05
tv, yv = integrate_sStepControl(f=f_1, t0=0.0, y0=Y0_f1, tend=100.0, h=h, tol=1.0E-05)
print("[ %20.15f, %20.15f]"%(yv[-1,0], yv[-1,1]) )
plt.plot(tv, yv)
得到上面的错误,但它被绘制出来了。不知道我在这里做错了什么:-/
编辑: 添加了对 err == 0
的检查您需要实际传回计算出的新步长 h
以便它可以在第一步的主循环中使用。
误差应按标准计算。添加一些小数以避免被零除。
四阶方法的主要误差项是 C*h^5
。这必须与所需的本地错误 tol*h
进行比较。总的来说,这会导致计算最优 h
的第 4 个根。取 5 次方根可以提供某种抑制作用,但是对全局误差的影响并不是很直接。
def rkck(f, x, y, h, tol):
#xn = x + h
err = 2 * tol
while (err > tol):
xn = x + h
k1 = h*f(x,y)
k2 = h*f(x+(1/5)*h,y+((1/5)*k1))
k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
dy4 = ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
dy5 = ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
err = 1e-2*tol+max(abs(dy4-dy5))
# h = 0.95 * h * (tol/err)**(1/5)
h = 0.8 * h * (tol*h/err)**(1/4)
yn = y+dy4
return xn, yn, h
def integrate_sStepControl(f, t0, y0, tend, h, tol):
T = [t0]
Y = [y0]
t = t0
y = y0
while (t < tend):
h = min(h, tend-t)
t, y, h = rkck(f, t, y, h, tol)
T.append(t)
Y.append(y)
return np.array(T), np.array(Y)
因此,您的示例给出了解决方案、时间步长和 error/tol
的以下图