Траектории фазового пространства для модели Хиндмарша-Роуза

#python #matplotlib #plot #numerical-methods #integrator

Вопрос:

Я пытаюсь построить траектории фазового пространства для модели Хиндмарша-Роуза. Я реализовал интегратор RK4 для решения следующего набора уравнений:

введите описание изображения здесь

Код, который я написал до сих пор, приведен ниже.

 import numpy as np
import matplotlib.pyplot as plt

def RK4(f, x0, t):
    
    dt   = t[2] -t[1]  #time span 
    N    = len(t)
    X    = np.empty((len(t), len(x0)))
    X[0] = x0
    
    for i in range(1, N):
        
        k1 = f(X[i-1], t[i-1])
        k2 = f(X[i-1]   dt/2*k1, t[i-1]   dt/2)
        k3 = f(X[i-1]   dt/2*k2, t[i-1]   dt/2)
        k4 = f(X[i-1]   dt*k3, t[i-1]   dt)
        
        X[i] = X[i-1]   dt/6*(k1   2*k2   2*k3   k4)
        
    return X


def hindmarsh(X, t):
    
    a  = 3.0
    c  = 1.0
    d  = 5.0
    s  = 4.0
    x0 = - 1.6
    
    # Bifurcation parameters
    
    b   = 3.09
    I   = 3.2
    eps = 0.001
      
    x,y,z = X
    
    dxdt  = y - (a * x**3)   (b * x**2)   I - z 
    dydt  = c - (d * x**2) - y 
    dzdt  = eps * ( (s * (x - x0)) - z) 
    
    return np.array([dxdt, dydt, dzdt])


T     = np.linspace(0,100,10000)

Y     = [0.03, 0.03, 3]

param = RK4( hindmarsh, Y, T )

ax = plt.axes(projection='3d')

zline = param[2]
yline = param[1]
xline = param[0]


ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

ax.plot3D(xline, yline, zline)

 

Однако вместо того, чтобы получить орбиту в фазовом пространстве, как показано на рисунке ниже, я получаю прямую линию через фазовое пространство. Я был бы признателен за любые советы о том, как получить сюжет ниже.

введите описание изображения здесь

Комментарии:

1. z является ли последний компонент, почему вы переключаете переменные при извлечении компонентов из решения? Постоянная Липшица составляет около 10, поэтому размер шага 0,01 находится на верхней границе допустимой точности.

2. @LutzLehmann это была ошибка, которую я допустил, когда пытался воспроизвести график. Спасибо, что заметили это. Я внес изменения в код.

Ответ №1:

param имеет форму (len(T), len(Y)) , поэтому время находится в первом измерении,а x,y, z-во втором измерении. Вы получите правильный сюжет с

 zline = param[:,0]
xline = param[:,1]
yline = param[:,2]