用Runge-Kutta算法求解Roessler振子

Question

我在这个问题上遇到了一些麻烦。我得到了这组方程

dx / dt = -y -z
dy / dt = x + a * y
dz / dt = b + z * (x - c)

以及默认值a=0.1、b=0.1、c=14以及龙格-库塔算法：

def rk4(f, xvinit, Tmax, N):
        T = np.linspace(0,Tmax,N+1)
        xv = np.zeros( (len(T), len(xvinit)) )
        xv[0] = xvinit
        h = Tmax / N
        for i in range(N):
            k1 = f(xv[i])
            k2 = f(xv[i] + h/2.0*k1)
            k3 = f(xv[i] + h/2.0*k2)
            k4 = f(xv[i] + h*k3)
            xv[i+1] = xv[i] + h/6.0 *( k1 + 2*k2 + 2*k3 + k4)
        return T, xv

我需要用0.1的时间步长求解这个从t=0到t=100的系统，使用初始条件(0,0,0)=(0,0,0) at =0我不太确定从哪里开始，我试着定义一个函数来给出振荡器：

def roessler(xyx, a=0.1, b=0.1, c=14):
    xyx=(x,y,x)
    dxdt=-y-z
    dydt=x+a*y
    dzdt=b+z*(x-c)
    return dxdt ,dydt ,dzdt       

它返回方程右侧的默认值，然后我尝试用roessler替换f，并用我给出的值填充xvinit，Tmax和N的值，但它不起作用。感谢任何帮助，如果这部分格式有误，我很抱歉，我是新来的。

Answer 1

好吧，你差不多已经明白了。将您的roessler函数更改为以下内容

def roessler(xyx, a=0.1, b=0.1, c=14):
    x, y, z = xyx
    dxdt=-y-z
    dydt=x+a*y
    dzdt=b+z*(x-c)
    return np.array([dxdt, dydt, dzdt])

然后调用

T, sol = rk4(roessler, np.array([0, 0, 0]), 100, 1000)

让它工作起来。

抛开roessler函数第一行中的拼写错误不谈，解决这个问题的关键是要理解您有一个微分方程的system，即，您需要使用vectors。虽然您已经将输入作为正确的向量，但您还需要使roessler的输出成为一个向量，并输入具有适当形状的初始值。