# Drawing saddle node bifurcation diagram for a non-linear ODE in Python

I'm trying to draw the bifurcation diagram of the following ODE, However what I get is not exactly right. There's a lot of "noise" as you can see on the figure below. Normally there should be the blue line (stable line) that goes from bottow left unto intersection point with orange line. Then the orange line that goes from blue line to green line (the one in the middle -- unstable line). Then the green line that goes from the intersection point with the orange line unto top right (unstable line).

In fact, it should be more like this (except it's the other way around but you can see the idea).

Here are my questions :

1. Any idea to fix / improve my algorithm ?
2. Is there a way I could precisely indicate the stable / unstable lines in matplotlib (e.g. red for unstable, green for stable) instead of all these colors ?
3. Is there any way to draw the flow of the dynamical system between the lines ?

Here is the code :

The idea of this code is

• Loop over all forcing paramaters (phi) of a forcing paramter range.
• For each forcing parameter, looking for some roots of the ODE.
• Plotting.
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 ** 3) + a2 * v + 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)

np.array([equilibria])
#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()



For the first question:

You should look into continuation. I don't know python so i can't write the code for you, but basically, you change the algorithm into something like this:

1. search for some solutions for phi = -2.
2. use those solutions as initial guesses when looking for phi = -1.9
3. repeat step 2 for the rest of the interval

Because you are using the previous solutions as guesses you will avoid branch jumping (the noise you spoke of). There are a lot of issues to take care of though. Continuation is a difficult algorithm to make work. I suggest you look for a library if you run into more trouble.

For the second question:

you can look at the derivative of $$\frac{dx}{dt}$$. If it is positive -> unstable, negative -> stable. Then use that information to draw the point in the correct color. I recommend using a dotted plot to facilitate this, looks better that way IMO.

For the third question:

No idea

• For (1), switching to a problem of maximisation and use particle swarm optimization algorithm (PSO) would be a good idea ? Seems like a fairly simple algorithm that does not involve much errors. Mar 3, 2021 at 8:17
• not really familiar with PSO, but sure. Although you'll lose the information about branching just like in your previous attempt. Mar 3, 2021 at 13:43