计算科学 - 在 Python 中绘制非线性 ODE 的鞍节点分岔图 - 吾爱随笔录

我正在尝试绘制以下 ODE 的分岔图，

此 ODE 导致鞍节点分叉 （参见 wiki）

但是我得到的并不完全正确。如下图所示，有很多“噪音”。

通常应该有从左下角到与橙色线交点的蓝线（稳定线）。然后是从蓝线到绿线的橙色线（中间的那条线——不稳定线）。然后是从与橙色线的交点到右上角的绿线（不稳定线）。

事实上，它应该更像这样（除了它是相反的，但你可以看到这个想法）。

这是我的问题：

有什么想法可以修复/改进我的算法吗？
有没有办法可以精确地指示 matplotlib 中的稳定/不稳定线（例如红色表示不稳定，绿色表示稳定）而不是所有这些颜色？
有什么方法可以在线条之间绘制动力系统的流动？

这是代码：

这段代码的想法是

循环遍历强制参数范围的所有强制参数 (phi)。
对于每个强制参数，寻找 ODE 的一些根。
绘图。

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import fsolve

# 1D EDO for Saddle-Node bifurcation paramaters
a1 = -1
a2 = 1

# initial conditions
x0 = 0.0
y0 = 5.0
z0 = 5.0

r0 = 5.0

# time range
t_init = 0
t_fin = 500
time_step = 0.01

def fold(v, phi):
  return np.array([
    a1 * (v[0] ** 3) + a2 * v[0] + phi
  ])


phi = 0
nphi = 100
nguesses = 3
phi_mesh = np.linspace(start=-2, stop=2, num=nphi)

# number of time steps
nt = int((t_fin - t_init) / time_step)
time_mesh = np.linspace(start=t_init, stop=t_fin, num=nt)

def fold_bifurcation():
    equilibria_mesh = np.zeros((nphi, nguesses))

    # for each phi
    for phi_index in range(0, nphi-1):

      # find the equilibria of our system
      guesses = np.linspace(start=-3, stop=3, num=nguesses)
      equilibria = []
      
      # look for some equilibria
      for guess in guesses:  
          equilibrium = fsolve(func=fold, x0=[guess], args=(phi_mesh[phi_index]))
          equilibria.append(equilibrium[0])
      
      np.array([equilibria])
      # add to the mesh
      #print(np.shape(equilibria_mesh[phi_index]), np.shape(equilibria))
        
      equilibria_mesh[phi_index] = np.array([equilibria])

    plot(dataset=equilibria_mesh.copy(), ylabel="x")

def plot(dataset, ylabel):
    plt.plot(phi_mesh, dataset)
    plt.xlabel("$\phi$")
    plt.ylabel(ylabel)
    plt.xlim(-2,2)
    plt.ylim(-5,5)
    plt.show()
    
fold_bifurcation()