# Intermediate tutorial: Computing stable and unstable branches of a 2D saddle point

(Note as of March 2015: This script only works with the latest git version of PyDSTool.)

A second order system governing the motion of a bead in a rotating hoop is

\left\{ \begin{align*} \frac{d\phi}{dt} &= \nu, \\ \frac{d\nu}{dt} &= -d \nu + g \sin(\phi)*\cos(\phi) - \sin(\phi) \end{align*} \right. \qquad (1)

where d and g are parameters, and $\phi$ is the phase, i.e. angular position.

This system is specified using PyDSTool with


import PyDSTool as dst
from PyDSTool.Toolbox import phaseplane as pp
import numpy as np
from matplotlib import pyplot as plt

# user selection
print("Trying to use Radau integrator, but use Vode")print("  if you don't have the C integrators working")

def build_sys():
# we must give a name

# parameters
DSargs.pars = {'g': 0,
'd': 0.3}

# rhs of the differential equation
DSargs.varspecs = {'phi': 'nu',
'nu': '-d*nu + g*sin(phi)*cos(phi) - sin(phi)'}

# initial conditions
DSargs.ics = {'phi': 0, 'nu': 0}

# set the domain of integration.
# (increased domain size to explore around phi=-pi saddle)
DSargs.xdomain = {'phi': [-2*np.pi, 2*np.pi], 'nu': [-4, 4]}

# allow tdomain to be infinite, set default tdata here
DSargs.tdata = [0, 50]

# to avoid typos / bugs, use built-in Symbolic differentation!
f = [DSargs.varspecs['phi'], DSargs.varspecs['nu']]
Df=dst.Diff(f, ['phi', 'nu'])
DSargs.fnspecs = {'Jacobian': (['t','phi','nu'],
str(Df.renderForCode()))}
# yields """[[0, 1], [g*cos(phi)*cos(phi) - g*sin(phi)*sin(phi) - cos(phi), -d]]""")}



### Auxiliary definitions

The system in Eq. (1) has multiple symmetries, and an elementary text on dynamical systems will explain why there are two saddle points and a stable spiral fixed point in the standard domain for $\phi \in \left[-\pi, \pi\right]$. We will not repeat that analysis here, and we will henceforth assume you know how to identify the non-degenerate, hyperbolic saddle point by its eigenvalues, and use its eigenvectors to locally approximate the stable and unstable sub-manifolds around it. (See this article for some more details and references.)

For PyDSTool to compute these manifolds, we also have to declare some auxiliary functions and events, which continue the build_sys function...


# Make auxiliary functions to define event lines near saddle
res = pp.make_distance_to_line_auxfn('Gamma_out_plus',
'Gamma_out_plus_fn',
('phi','nu'), True)
man_pars = res['pars']
man_auxfns = res['auxfn']
res = pp.make_distance_to_line_auxfn('Gamma_out_minus',
'Gamma_out_minus_fn',
('phi','nu'), True)
man_pars.extend(res['pars'])
man_auxfns.update(res['auxfn'])

# update param values with defaults (0)
for p in man_pars:
DSargs.pars[p] = 0

if gentype in [dst.Generator.Vode_ODEsystem, dst.Generator.Euler_ODEsystem]:
targetlang = 'python'
else:
targetlang = 'c'

DSargs.fnspecs.update(man_auxfns)
ev_plus = dst.Events.makeZeroCrossEvent(expr='Gamma_out_plus_fn(%s,%s)'%('phi','nu'),
dircode=0,
argDict={'name': 'Gamma_out_plus',
'eventtol': 1e-5,
'eventdelay': 1e-3,
'starttime': 0,
'precise': False,
'active': False,
'term': True},
targetlang=targetlang,
varnames=['phi','nu'],
fnspecs=man_auxfns,
parnames=man_pars
)
ev_minus = dst.Events.makeZeroCrossEvent(expr='Gamma_out_minus_fn(%s,%s)'%('phi','nu'),
dircode=0,
argDict={'name': 'Gamma_out_minus',
'eventtol': 1e-5,
'eventdelay': 1e-3,
'starttime': 0,
'precise': False,
'active': False,
'term': True},
targetlang=targetlang,
varnames=['phi','nu'],
fnspecs=man_auxfns,
parnames=man_pars
)

DSargs.events = [ev_plus, ev_minus]

# an instance of the 'Generator' class.
return gentype(DSargs)




PyDSTool uses a naive algorithm for finding the sub-manifolds. It is a "quick and dirty" approach that could easily be optimized for specific systems and situations. It may also not work well for very stiff systems, as the backwards integration necessary to find the unstable sub-manifolds may not be feasible.

### Basic phase plane information

The tutorial on the van der Pol oscillator showed the basics of displaying phase plane information, which we will not repeat here. This is the code to bring up the basic information, and we don't repeat the custom phaseplane plotting functions mentioned here (see the actual code for details).


# plot vector field, using a scale exponent to ensure arrows are well spaced
# and sized
plot_PP_vf_custom(ode_sys, 'phi', 'nu', scale_exp=-0.25)

# find fixed points
fp_coords = pp.find_fixedpoints(ode_sys, eps=1e-6)

# n=3 uses three starting points in the domain to find nullcline parts, to an
# accuracy of eps=1e-8, and a maximum step for the solver of 0.1 units.
# The fixed points found is also provided to help locate the nullclines.
nulls_x, nulls_y = pp.find_nullclines(ode_sys, 'phi', 'nu', n=3,
eps=1e-6, max_step=0.1, fps=fp_coords)

# plot the fixed points
fps = []
for fp_coord in fp_coords:
fps.append( pp.fixedpoint_2D(ode_sys, dst.Point(fp_coord)) )

for fp_obj in fps:
plot_PP_fps_custom(fp_obj, do_evecs=True, markersize=7, flip_coords=True)

# plot the nullclines
plt.plot(nulls_x[:,0], nulls_x[:,1], 'b')
plt.plot(nulls_y[:,0], nulls_y[:,1], 'g')

plt.axis('tight')
plt.title('Phase plane')
plt.xlabel('phi')
plt.ylabel('nu')

# you may not need to run these commands on your system
plt.draw()
plt.show()



We set up some internal parameters for the integrator, then define the saddle point (we already know it is index 1 of the list of fixed points, fps, and the function find_saddle_manifolds will verify its classification too). We build the manifold parts up as a dictionary keyed by the stable 's' or unstable 'u' characters. Each of those contains a dictionary keyed by direction (arbitrary which is which) given by +/- 1 (integer). For all four, we begin by making a singleton Pointset object containing the saddle point itself. Those Pointsets will be parameterized by arclength relative to zero at the saddle point. This makes positive arclength in the 'forwards' direction away from the saddle, and negative arclengths in the 'backwards' direction.


# magBound change ensures quicker determination of divergence during
# manifold computations. max_pts must be larger when we are further
# away from the fixed point.
ode_sys.set(algparams={'magBound': 10000})

# KEY to manifold parts dictionary key codes
# u means unstable branch
# s means stable branch
# +1 means 'forwards' direction
# -1 means 'backwards' direction

manifold_parts = {
'u': {1: dst.Pointset(indepvarname='arc_len',
indepvararray=[0],
-1: dst.Pointset(indepvarname='arc_len',
indepvararray=[0],
's': {1: dst.Pointset(indepvarname='arc_len',
indepvararray=[0],
-1: dst.Pointset(indepvarname='arc_len',
indepvararray=[0],
}

# t = 60 is probably overkill but the events are terminal so won't compute more than needed


Next, we compute each branch in two stages. We can afford much smaller initial perturbations close to the saddle point, which will make the computation of the first few points much quicker in the first stage. The second stage must use a wider range to correct the initial guess because the manifold is probably more curved and divergent from the linear local sub-manifold. The two variables are of comparable scale so we keep the default argument rel_scale as 1:1. The other parameters are described graphically in the figure below.

In the schematic diagram, the saddle point is in the center as an open circle. The thick green (red) arrows represent the eigenvectors in the stable (unstable) direction. These are the local approximations to the sub-manifolds, $W^s_{loc}$ and $W^u_{loc}$. Tangent to these at the saddle point are the actual sub-manifolds, $W^s$ and $W^u$, shown by the thinner line of the same color. These are typically curved. The diagram summarizes the initial and next step of the calculation for one branch of the stable sub-manifold. The two event thresholds are given by $\Gamma_{+}$ and $\Gamma_{-}$, and shown as thick black lines. These events are automatically generated based on the vector that's normal to the unstable direction's eigenvector, at a distance set by the ds_gamma parameter.

Two example forward trajectories are shown in blue, diverging either side of the saddle point. In this case, each one crosses the event thresholds at points shown by the yellow stars.

The ic_ds parameter is the offset from the saddle point along the stable eigenvector where we start the predictor-corrector process. The initial guess (result of the prediction step) is then at the square-marked point $\tilde{x}_0$. It is corrected by a simple bisection step, given that we assume the solution is bounded above and below by the points $a=\tilde{x}_0+\mathrm{ds\_perp} . \hat{n}, \ b=\tilde{x}_0-\mathrm{ds\_perp} . \hat{n}$, where $\hat{n}$ is the normal direction to the eigenvector. Assuming we choose ds_perp large enough to ensure that the true solution $x_0$ lies in $[a,b]$, then bisection will find it to a required tolerance (eps). For the iterative step that creates the next predicted point on the manifold, we move an amount of arclength ds backwards along the flow from $x_0$. We then repeat the correction steps as above.

Note that, if ds_perp is too large, then we may fail to cross one of the event thresholds before we run out of integrated points. If that happens, ds_perp will be automatically reduced in size (default by a factor of 0.75) until one of the targets is reached. If that doesn't happen before the tolerance ds_perp_eps is reached, then the algorithm will fail.

The code to compute the manifolds looks like this:


verbose = 0
max_pts = 600
for which_man in ['s', 'u']:
for dirn in [-1, 1]:
man_part = manifold_parts[which_man][dirn]
print("Starting manifold %s, direction %i"%(which_man, dirn))
ds_perp = 0.02 # initial value for stage 2
ds_gamma = 0.3 # initial value for stage 2
# for speed, and because of symmetry, only compute long arcs on one side
if dirn == 1:
max_arclen = 4
else:
max_arclen = 9

while len(man_part) < max_pts and max(abs(man_part['arc_len'])) < max_arclen:
attempt_num = 0
# Perform calculation in two stages. First stage is more accurate while we
# are close to the saddle point
if len(man_part) == 1:
# first stage (only called once)
print "  First stage..."
ode_sys.set(algparams={'max_pts': 20000})
ds_perp=0.005, tmax=60, max_arclen=max_arclen, eps=2e-5,
ic_ds=0.0002, max_pts=250, directions=(dirn,), ev_dirn=1,
which=(which_man,), other_pts=[fps[0].point, fps[2].point],
rel_scale=(1,1), verboselevel=verbose, fignum=1)
part = man_new[which_man][dirn]
if dirn == 1:
select = -1
else:
# backwards means that points will be added behind the initial point
# in terms of Pointset order, so the "last" point computed will be
# at index 0
select = 0
part.indepvararray += man_part.indepvararray[select]
man_part.insert(part)
attempt_num = 1
else:
# Stage two (called repeatedly until length while loop satisfied)
print("  Continuing in stage 2...")
ode_sys.set(algparams={'max_pts': 100000})
if dirn == 1:
select = -1
else:
# backwards means that points will be added behind the initial point
# in terms of Pointset order, so the "last" point computed will be
# at index 0
select = 0
attempt_num += 1
try:
# groups of max_pts/4 at a time
ds_perp=ds_perp, tmax=40, max_arclen=max_arclen, eps=2e-5,
ic=man_part[select], max_pts=int(max_pts/4.0),
directions=(dirn,), ev_dirn=1,
which=(which_man,), other_pts=[fps[0].point, fps[2].point],
rel_scale=(1,1), verboselevel=verbose, fignum=1)
except RuntimeError:
# proceed with what we've got
print("Initial convergence error: Proceeding with what we've got!")
ds_perp *= 2
ds_gamma *= 4
break # to continue
else:
part = man_new[which_man][dirn]
part.indepvararray += man_part.indepvararray[select]
man_part.insert(part)

plot_manifold(manifold_parts, 's', 'k-')
plot_manifold(manifold_parts, 'u', 'r-')


Unfortunately, we have to care about the directions that we are allowed to cross $\Gamma_{\pm}$ because, despite the diagram's suggestion of their finite length, their representation as a zero-crossing event hypersurface (line, in this case) in the code means they extend infinitely in both directions. So, when we are far from the saddle our starting points may cross these surfaces in the 'wrong' direction. To realize which direction we have to set the event detection parameter ev_dirn, we can either use trial and error or figure it out by thinking about the orientation of the eigenvectors and whether we are studying the stable manifold (in forward time) or the unstable manifold (in backwards time). See the function's docstring for more information.

This tutorial code yields the following diagram. The thin blue and green lines are the nullclines. The red and black lines are the unstable and stable sub-manifolds of the saddle, respectively. You can see how they follow the flow lines shown by the quiver plot.

Thanks to Brian Merchant for some of the materials used to make this tutorial.