如何阻止我的 Runge-Kutta2 (Heun) 方法爆炸?
How can I stop my Runge-Kutta2 (Heun) method from exploding?
我目前正在尝试编写一些 python 代码来求解一阶 ODE 的任意系统,使用由值 alpha、gamma(两个维度 m 的向量)定义的一般显式 Runge-Kutta 方法和用户传入的 Butcher table 的 beta(维度为 m x m 的下三角矩阵)。我的代码似乎适用于单个 ODE,已经在几个不同的示例中对其进行了测试,但我正在努力将我的代码推广到向量值 ODE(即系统)。
特别是,我尝试使用代码中给出的 Butcher Tableau 值定义的 Heun 方法求解 Van der Pol 振荡器 ODE(简化为一阶系统),但我收到错误
- "RuntimeWarning: overflow encountered in double_scalars
f = lambda t,u: np.array(... etc)
" and
- "RuntimeWarning: invalid value encountered in add
kvec[i] = f(t+alpha[i]*h,y+h*sum)
"
其次是我的解决方案向量,它显然已经爆炸了。请注意,下面注释掉的代码是我尝试并正确求解的单个 ODE 示例之一。有人可以帮忙吗?这是我的代码:
import numpy as np
def rk(t,y,h,f,alpha,beta,gamma):
'''Runga Kutta iteration'''
return y + h*phi(t,y,h,f,alpha,beta,gamma)
def phi(t,y,h,f,alpha,beta,gamma):
'''Phi function for the Runga Kutta iteration'''
m = len(alpha)
count = np.zeros(len(f(t,y)))
kvec = k(t,y,h,f,alpha,beta,gamma)
for i in range(1,m+1):
count = count + gamma[i-1]*kvec[i-1]
return count
def k(t,y,h,f,alpha,beta,gamma):
'''returning a vector containing each step k_{i} in the m step Runga Kutta method'''
m = len(alpha)
kvec = np.zeros((m,len(f(t,y))))
kvec[0] = f(t,y)
for i in range(1,m):
sum = np.zeros(len(f(t,y)))
for l in range(1,i+1):
sum = sum + beta[i][l-1]*kvec[l-1]
kvec[i] = f(t+alpha[i]*h,y+h*sum)
return kvec
def timeLoop(y0,N,f,alpha,beta,gamma,h,rk):
'''function that loops through time using the RK method'''
t = np.zeros([N+1])
y = np.zeros([N+1,len(y0)])
y[0] = y0
t[0] = 0
for i in range(1,N+1):
y[i] = rk(t[i-1],y[i-1], h, f,alpha,beta,gamma)
t[i] = t[i-1]+h
return t,y
#################################################################
'''f = lambda t,y: (c-y)**2
Y = lambda t: np.array([(1+t*c*(c-1))/(1+t*(c-1))])
h0 = 1
c = 1.5
T = 10
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])
eff_rk = compute(h0,Y(0),T,f,alpha,beta,gamma,rk, Y,11)'''
#constants
mu = 100
T = 1000
h = 0.01
N = int(T/h)
#initial conditions
y0 = 0.02
d0 = 0
init = np.array([y0,d0])
#Butcher Tableau for Heun's method
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])
#rhs of the ode system
f = lambda t,u: np.array([u[1],mu*(1-u[0]**2)*u[1]-u[0]])
#solving the system
time, sol = timeLoop(init,N,f,alpha,beta,gamma,h,rk)
print(sol)
您的步长不够小。 mu=100
的 Van der Pol 振荡器是一个快慢系统,在模式切换时有非常急转弯,因此相当僵硬。对于显式方法,这需要较小的步长,最小的合理步长是 1e-5
到 1e-6
。您已经得到了 h=0.001
的极限环解,结果速度高达 150。
您可以通过使用不同的 velocity/impulse 变量来降低一些刚度。在等式中
x'' - mu*(1-x^2)*x' + x = 0
你可以把前两项组合成导数,
mu*v = x' - mu*(1-x^2/3)*x
这样
x' = mu*(v+(1-x^2/3)*x)
v' = -x/mu
第二个方程现在在接近极限环时均匀变慢,而第一个方程在 v
离开立方体 v=x^3/3-x
时有较长的相对直线跳跃。
这与原始 h=0.01
很好地集成,将解决方案保持在框内 [-3,3]x[-2,2]
,即使它显示了一些奇怪的振荡,这些振荡对于较小的步长和精确的解决方案是不存在的。
我目前正在尝试编写一些 python 代码来求解一阶 ODE 的任意系统,使用由值 alpha、gamma(两个维度 m 的向量)定义的一般显式 Runge-Kutta 方法和用户传入的 Butcher table 的 beta(维度为 m x m 的下三角矩阵)。我的代码似乎适用于单个 ODE,已经在几个不同的示例中对其进行了测试,但我正在努力将我的代码推广到向量值 ODE(即系统)。
特别是,我尝试使用代码中给出的 Butcher Tableau 值定义的 Heun 方法求解 Van der Pol 振荡器 ODE(简化为一阶系统),但我收到错误
- "RuntimeWarning: overflow encountered in double_scalars
f = lambda t,u: np.array(... etc)
" and- "RuntimeWarning: invalid value encountered in add
kvec[i] = f(t+alpha[i]*h,y+h*sum)
"
其次是我的解决方案向量,它显然已经爆炸了。请注意,下面注释掉的代码是我尝试并正确求解的单个 ODE 示例之一。有人可以帮忙吗?这是我的代码:
import numpy as np
def rk(t,y,h,f,alpha,beta,gamma):
'''Runga Kutta iteration'''
return y + h*phi(t,y,h,f,alpha,beta,gamma)
def phi(t,y,h,f,alpha,beta,gamma):
'''Phi function for the Runga Kutta iteration'''
m = len(alpha)
count = np.zeros(len(f(t,y)))
kvec = k(t,y,h,f,alpha,beta,gamma)
for i in range(1,m+1):
count = count + gamma[i-1]*kvec[i-1]
return count
def k(t,y,h,f,alpha,beta,gamma):
'''returning a vector containing each step k_{i} in the m step Runga Kutta method'''
m = len(alpha)
kvec = np.zeros((m,len(f(t,y))))
kvec[0] = f(t,y)
for i in range(1,m):
sum = np.zeros(len(f(t,y)))
for l in range(1,i+1):
sum = sum + beta[i][l-1]*kvec[l-1]
kvec[i] = f(t+alpha[i]*h,y+h*sum)
return kvec
def timeLoop(y0,N,f,alpha,beta,gamma,h,rk):
'''function that loops through time using the RK method'''
t = np.zeros([N+1])
y = np.zeros([N+1,len(y0)])
y[0] = y0
t[0] = 0
for i in range(1,N+1):
y[i] = rk(t[i-1],y[i-1], h, f,alpha,beta,gamma)
t[i] = t[i-1]+h
return t,y
#################################################################
'''f = lambda t,y: (c-y)**2
Y = lambda t: np.array([(1+t*c*(c-1))/(1+t*(c-1))])
h0 = 1
c = 1.5
T = 10
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])
eff_rk = compute(h0,Y(0),T,f,alpha,beta,gamma,rk, Y,11)'''
#constants
mu = 100
T = 1000
h = 0.01
N = int(T/h)
#initial conditions
y0 = 0.02
d0 = 0
init = np.array([y0,d0])
#Butcher Tableau for Heun's method
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])
#rhs of the ode system
f = lambda t,u: np.array([u[1],mu*(1-u[0]**2)*u[1]-u[0]])
#solving the system
time, sol = timeLoop(init,N,f,alpha,beta,gamma,h,rk)
print(sol)
您的步长不够小。 mu=100
的 Van der Pol 振荡器是一个快慢系统,在模式切换时有非常急转弯,因此相当僵硬。对于显式方法,这需要较小的步长,最小的合理步长是 1e-5
到 1e-6
。您已经得到了 h=0.001
的极限环解,结果速度高达 150。
您可以通过使用不同的 velocity/impulse 变量来降低一些刚度。在等式中
x'' - mu*(1-x^2)*x' + x = 0
你可以把前两项组合成导数,
mu*v = x' - mu*(1-x^2/3)*x
这样
x' = mu*(v+(1-x^2/3)*x)
v' = -x/mu
第二个方程现在在接近极限环时均匀变慢,而第一个方程在 v
离开立方体 v=x^3/3-x
时有较长的相对直线跳跃。
这与原始 h=0.01
很好地集成,将解决方案保持在框内 [-3,3]x[-2,2]
,即使它显示了一些奇怪的振荡,这些振荡对于较小的步长和精确的解决方案是不存在的。