## Saddle node bifurcation at infinity

Author: GC

Date: Dec. 10, 2016

Geographical data visualization

The system defined by $\begin{cases} \dot{x} = 1 + \varepsilon x^2 \\ \dot{y} = -y \end{cases}$ admits a saddle node bifurcation at infinity. The next program produces a nice view of that bifurcation on the Poincaré sphere.

### Source code

#!/usr/bin/env python3

"""
Program to visualize a saddle node bifurcation at infinity on the Poincaré sphere
"""

# Scientific libraries
from matplotlib import pyplot as plt
from matplotlib.pyplot import cm
import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
import matplotlib.gridspec as gridspec
from matplotlib.patches import Ellipse, PathPatch, Circle
import mpl_toolkits.mplot3d.art3d as art3d

cmap = plt.get_cmap('Blues')
def trajectory(IV, until, t0):
orbit = odeint(system, IV, np.arange(0, until, 0.01))
x, y = orbit.T
X = x/np.sqrt(1+x**2+y**2)
Y = y/np.sqrt(1+x**2+y**2)
Z = 1/np.sqrt(1+x**2+y**2)
ax.plot3D(X, Y, Z, color = cmap(0.7-(+X[0]+X[-1])/5), linewidth=1.5)
a = Arrow3D([X[t0],X[t0+2]],[Y[t0],Y[t0+2]],[Z[t0],Z[t0+2]],
mutation_scale=20,
lw=1, arrowstyle="->",
color=cmap(0.7-(X[0]+X[-1])/5))

def sphere(cstride_wireframe):
ax.plot_surface(x1, y1, z1, rstride=1, cstride=1, linewidth=0.0, color='#FAFAFA')

def equator():
ax.plot(x2, y2, z2, lw=2.5, color='brown', alpha=0.7)

def view(e, a):
ax.imshow([[0, -1], [0, 1]],
interpolation='bicubic',
cmap=cm.winter,
extent=(-0.1, 0.1, -0.1, 0.1),
alpha=0.6)
ax.set_aspect('auto')
ax.view_init(elev=e, azim=a)

def limits(u):
ax.set_xlim(-u, u)
ax.set_ylim(-u, u)
ax.set_zlim(-u, u)

def equilibrium(eps):
x = np.sqrt(-1/eps)
y = 0
X = x/np.sqrt(1+x**2+y**2)
Y = y/np.sqrt(1+x**2+y**2)
Z = 1/np.sqrt(1+x**2+y**2)
point1 = Ellipse((X, Y), 0.05, 0.06, angle=25)
Ellipse.set_color(point1,'0.2')
art3d.pathpatch_2d_to_3d(point1, z=Z, zdir="z")
x = -np.sqrt(-1/eps)
y = 0
X = x/np.sqrt(1+x**2+y**2)
Y = y/np.sqrt(1+x**2+y**2)
Z = 1/np.sqrt(1+x**2+y**2)
point2 = Ellipse((X, Y), 0.1, 0.07, angle=-40)
Ellipse.set_color(point2,'white')
art3d.pathpatch_2d_to_3d(point2, z=Z, zdir="z")

# 3d arrow
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs

def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)

# Sphere
u = np.linspace(0, 2 * np.pi, 39)
v = np.linspace(0, np.pi, 21)
x1 = np.outer(np.cos(u), np.sin(v))
y1 = np.outer(np.sin(u), np.sin(v))
z1 = np.outer(np.ones(np.size(u)), np.cos(v))

# Equator
u = np.linspace(0, 2 * np.pi, 100)
x2 = np.cos(u)
y2 = np.sin(u)
z2 = np.array([0 for k in range(len(u))])

#System
alpha = 1

def system(X, t):
x, y = X
dx = alpha + eps*x**2
dy = -y
return [dx, dy]

# Figure
plt.figure(figsize=(10,10))
gs = gridspec.GridSpec(2, 2, wspace=0.0, hspace=0.0)

# Poincaré sphere 1
eps = -0.5
ax = plt.subplot(gs[0, 0], projection='3d')
ax.set_axis_off()
ax.grid(False)
limits(0.6)
sphere(4)
equator()
trajectory([-1, 100], 6.5, 350)
trajectory([-1, -100], 6.5, 350)
trajectory([100, 0], 2, 100)
trajectory([-1.61, 0], 1.92, 130)
trajectory([-1.3, 0], 4.3, 250)
trajectory([-np.sqrt(2), 100], 6.5, 400)
trajectory([-np.sqrt(2), -100], 6.5, 400)
trajectory([-1.4, 100], 6.5, 450)
trajectory([-1.4, -100], 6.5, 450)
trajectory([-1.41, 100], 6.5, 450)
trajectory([-1.41, -100], 6.5, 450)
trajectory([-1.4125, 100], 6.5, 450)
trajectory([-1.4125, -100], 6.5, 450)
trajectory([-1.415, 100], 5.75, 450)
trajectory([-1.415, -100], 5.75, 450)
trajectory([-1.4138, 100], 7, 550)
trajectory([-1.4138, -100], 7, 550)
trajectory([-1.43, 100], 3.67, 320)
trajectory([-1.43, -100], 3.67, 320)
trajectory([200, 5], 3, 150)
trajectory([200, -5], 3, 150)
equilibrium(eps)
view(-21, -36)

# Poincaré sphere 2
eps = -0.1
ax = plt.subplot(gs[0, 1], projection='3d')
ax.set_axis_off()
ax.grid(False)
limits(0.6)
sphere(4)
equator()
trajectory([-1, 100], 6.5, 350)
trajectory([-1, -100], 6.5, 350)
trajectory([-2.6, 0], 6.5, 200)
trajectory([-4.2, 0], 2.8, 120)
trajectory([100, 0], 4, 200)
trajectory([-3.1, 100], 10.5, 520)
trajectory([-3.1, -100], 10.5, 520)
trajectory([-3, 100], 8.5, 520)
trajectory([-3, -100], 8.5, 520)
trajectory([-2.7, 100], 6.5, 420)
trajectory([-2.7, -100], 6.5, 420)
trajectory([-2.5, 100], 6.5, 320)
trajectory([-2.5, -100], 6.5, 320)
trajectory([-2.9, 100], 6.5, 400)
trajectory([-2.9, -100], 6.5, 400)
trajectory([-np.sqrt(-1/eps), 100], 6.2, 400)
trajectory([-np.sqrt(-1/eps), -100], 6.2, 400)
trajectory([-np.sqrt(-1/eps)-0.1, 100], 6, 400)
trajectory([-np.sqrt(-1/eps)-0.1, -100], 6, 400)
equilibrium(eps)
view(-31, -26)

# Poincaré sphere 3
eps = -0.01
ax = plt.subplot(gs[1, 0], projection='3d')
ax.set_axis_off()
ax.grid(False)
limits(0.6)
sphere(3)
equator()
trajectory([-1, 100], 8, 350)
trajectory([-1, -100], 8, 350)
trajectory([-7, 0], 10, 700)
trajectory([-10, 100], 5.5, 300)
trajectory([-10, -100], 5.5, 300)
trajectory([-5, 100], 10, 400)
trajectory([-5, -100], 10, 400)
trajectory([-6, 100], 10, 400)
trajectory([-6, -100], 10, 400)
trajectory([-7, 100], 10, 400)
trajectory([-7, -100], 10, 400)
trajectory([-3, 100], 10, 400)
trajectory([-3, -100], 10, 400)
trajectory([-4, 100], 10, 400)
trajectory([-4, -100], 10, 400)
trajectory([10, 100], 5, 250)
trajectory([10, -100], 5, 250)
equilibrium(eps)
view(-41, -50)

# Poincaré sphere 4
eps=0.1
ax = plt.subplot(gs[1, 1], projection='3d')
ax.set_axis_off()
ax.grid(False)
limits(0.6)
sphere(3)
equator()
trajectory([-1, 100], 5.7, 350)
trajectory([-1, -100], 5.7, 350)
trajectory([-50, 0], 9.5, 450)
trajectory([-100, 20], 9.5, 450)
trajectory([-100, -20], 9.5, 450)
trajectory([-100, 60], 9.5, 450)
trajectory([-100, -60], 9.5, 450)
trajectory([-100, 100], 9.5, 450)
trajectory([-100, -100], 9.5, 450)
trajectory([-100, 200], 9.5, 450)
trajectory([-100, -200], 9.5, 450)
trajectory([-100, 400], 9.5, 550)
trajectory([-100, -400], 9.5, 550)
view(-51, -6)

plt.show()