如何阻止我的 Runge-Kutta2 (Heun) 方法爆炸? [英] How can I stop my Runge-Kutta2 (Heun) method from exploding?

问题描述

我目前正在尝试编写一些 python 代码来解决一阶 ODE 的任意系统，使用由值 alpha、gamma(两个维度为 m 的向量)和 beta(下三角矩阵)定义的通用显式 Runge-Kutta 方法由用户传入的 Butcher 表的维度 (mxm).我的代码似乎适用于单个 ODE，已经在几个不同的示例上对其进行了测试，但我正在努力将我的代码推广到向量值 ODE(即系统).

I am currently trying to write some python code to solve an arbitrary system of first order ODEs, using a general explicit Runge-Kutta method defined by the values alpha, gamma (both vectors of dimension m) and beta (lower triangular matrix of dimension m x m) of the Butcher table which are passed in by the user. My code appears to work for single ODEs, having tested it on a few different examples, but I'm struggling to generalise my code to vector valued ODEs (i.e. systems).

特别是，我尝试使用由我的代码中给出的 Butcher Tableau 值定义的 Heun 方法来解决 Van der Pol 振荡器 ODE(简化为一阶系统)，但我收到错误

In particular, I try to solve a Van der Pol oscillator ODE (reduced to a first order system) using Heun's method defined by the Butcher Tableau values given in my code, but I receive the errors

• "RuntimeWarning: 在 double_scalars 中遇到溢出 f = lambda t,u: np.array(... etc)" 和
• "RuntimeWarning: 在添加中遇到无效值 kvec[i] = f(t+alpha[i]*h,y+h*sum)"

• "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 示例之一.有人可以帮忙吗?这是我的代码:

followed by my solution vector that is clearly blowing up. Note that the commented out code below is one of the examples of single ODEs that I tried and is solved correctly. Could anyone please help? Here is my code:

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.

Your step size is not small enough. The Van der Pol oscillator with mu=100 is a fast-slow system with very sharp turns at the switching of the modes, so rather stiff. With explicit methods this requires small step sizes, the smallest sensible step size is 1e-5 to 1e-6. You get a solution on the limit cycle already for h=0.001, with resulting velocities up to 150.

您可以通过使用不同的速度/脉冲变量来降低一些刚度.在等式中

You can reduce some of that stiffness by using a different velocity/impulse variable. In the equation

    x'' - mu*(1-x^2)*x' + x = 0

你可以将前两项组合成一个导数，

you can combine the first two terms into a derivative,

    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 时，第一个方程有较长的相对直线跳跃代码>.

The second equation is now uniformly slow close to the limit cycle, while the first has long relatively straight jumps when v leaves the cubic v=x^3/3-x.

这与原始 h=0.01 很好地集成在一起，将解决方案保留在 [-3,3]x[-2,2] 框内，即使它显示了一些奇怪的振荡，这些振荡对于较小的步长和精确解是不存在的.

This integrates nicely with the original h=0.01, keeping the solution inside the box [-3,3]x[-2,2], even if it shows some strange oscillations that are not present for smaller step sizes and the exact solution.

这篇关于如何阻止我的 Runge-Kutta2 (Heun) 方法爆炸?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！