Участок бифуркации Хопфа

Question

Я пытаюсь закодировать бифуркационную диаграмму, чтобы проиллюстрировать значения f, для которых модель Oregonator дает колебательное поведение. Я получаю ошибку «установка элемента массива с последовательностью» в строке solve_ivp. Я подозреваю, что это имеет какое-то время, но я не уверен. Я должен получить бифуркацию Хопфа, то есть пулевидный конус в области колебаний.

Вот код:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Dimensionless constant parameters
eps = 0.04
a = 0.0008

# Dimensionless varying parameter - will reveal limit cycle region
f = np.linspace(-5,5,250)

# Oregonator model
def Oregonator(t, Y):
    x,z = Y;
    return [(x * (1 - x) + ((a - x) * f * z) / (a + x)) / eps, x - z]

# Time span, initial conditions

ts = np.linspace(-5, 5, 250)
Y0 = [1, 0.5]


# Numerical algorithm/method
NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
t = NumSol.t
x,z = NumSol.y

# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')


ax.plot(f, x, z, 'm')
ax.set_xlabel('$f$', fontsize=10)
ax.set_ylabel('$x^*$', fontsize=10)
ax.set_zlabel('$z^*$', fontsize=10)

ax.axis([-5, 5, -5, 5])
plt.grid()
plt.show() 

Answer 1

Вы не можете делать этот вид одновременного вычисления. (То есть, вы можете, но это требует явного кодирования и все еще не рекомендуется, поскольку выбор размера шага может сильно варьироваться в диапазоне f значений.)

Вы должны вычислить решение для каждого из f значения отдельно, а затем построить их. Чтобы сначала создать список всех решений, полезно собрать конструкцию для одного решения и его последней точки в отдельную функцию, чтобы построение списка можно было выполнить с помощью обработки списка.

# Dimensionless constant parameters
eps = 0.04
a = 0.0008

def limit(f):
    # Oregonator model
    def Oregonator(t, Y):
        x,z = Y;
        return [(x * (1 - x) + ((a - x) * f * z) / (a + x)) / eps, x - z]

    # Time span, initial conditions

    ts = np.linspace(-5, 5, 250)
    Y0 = [1, 0.5]


    # Numerical algorithm/method
    NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
    t = NumSol.t
    x,z = NumSol.y
    return x[-1],z[-1]

# Dimensionless varying parameter - will reveal limit cycle region
f = np.linspace(-5,5,250)

x,z = np.array([ limit(ff) for ff in f ]).T

# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot(f, x, z, 'mo', ms=2)

ax.set_xlabel('$f$', fontsize=10)
ax.set_ylabel('$x^*$', fontsize=10)
ax.set_zlabel('$z^*$', fontsize=10)

ax.axis([-5, 5, -5, 5])
plt.grid()
plt.show() 

Это результаты на участке