PyDSTool / Discussion / Help: Plotting stable and unstable manifolds of a 2D saddles?

Hi all,

I would like to plot the stable and unstable manifolds of 2D saddles.

In PyDSTool -> Toolbox -> phaseplane.py, I know there is the find_saddle_manifolds function:

def find_saddle_manifolds(fp, dx=None, dx_gamma=None, dx_perp=None, tmax=None,
                          max_len=None, ic=None,
                          ic_dx=None, max_pts=1000, directions=(1,-1),
                          which=('s', 'u'), other_pts=None, rel_scale=None,
                          dx_perp_fac=0.75, verboselevel=0, fignum=None):
    """Compute any branch of the stable or unstable sub-manifolds of a saddle.
    Accepts fixed point instances of class fixedpoint_2D.

    Required inputs:
      fp:       fixed point object
      dx:       arc-length step size (**fixed**)
      dx_gamma: determines the positions of the Gamma_plus and Gamma_minus
                event surfaces (can be a real scalar or a pair if not symmetric)
      dx_perp:  initial perturbation from the local linear sub-manifolds to
                find starting points.
      tmax:     maximum time to compute a trajectory before 'failing'
      max_len:  maximum arc length to compute
      max_pts:  maximum number of points to compute on each sub-manifold branch
      Specify either ic or ic_dx for initial point (e.g. to restart the calc
        after a previous failure) or a certain distance from the saddle point.

    Optional inputs:
      rel_scale:  a pair giving relative scalings of x and y coordinates in
        the plane, to improve stepping in the different directions.
        e.g. (1,10) would make dx steps in the y-direction 10 times larger than
        in the x-direction.

      which:  which sub-manifold to compute 's', 'u' or ('s', 'u').
        Default is both.

      directions:  which directions along chosen sub-manifolds? (1,), (-1,)
        or (1,-1). Default is both.

      other_pts can be a list of points whose proximity will be checked,
    and the computation halted if they get within dx of the manifold.

      dx_perp_fac:  For advanced use only. If you get failures saying dx_perp
         too small and that initial displacement did not straddle manifold, try
         increasing this factor towards 1 (default 0.75). Especially for
         unstable manifolds, initial values for dx_perp may diverge, but if
         dx_perp is shrunk too quickly with this factor the sweet spot may be
         missed.

      verboselevel
      fignum
    """

The problem is that I don't understand some of the required arguments:

dx_gamma: determines the positions of the Gamma_plus and Gamma_minus event surfaces (can be a real scalar or a pair if not symmetric)

dx_perp: initial perturbation from the local linear sub-manifolds to
find starting points

max_len: maximum arc length to compute

max_pts: maximum number of points to compute on each sub-manifold branch

Can you help me understand what these arguments mean, and what I should supply for them?

My second question is this: what sort of object is manifold (the object returned by the function)? How can I plot it?

For reference, here's my current PyDSTool code (plots phase plane, and fixed points of my dynamical system):

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

def plot_PP_fps_custom(fps, coords=None, do_evecs=False, markersize=10, flip_coords=False):
    """Draw 2D list of fixed points (singletons allowed), must be
    fixedpoint_2D objects.

    Optional do_evecs (default False) draws eigenvectors around each f.p.

    Requires matplotlib
    """
    if isinstance(fps, pp.fixedpoint_2D):
        # singleton passed
        fps = [fps]

    x, y = fps[0].fp_coords
    for fp in fps:
        # When matplotlib implements half-full circle markers
        #if fp.classification == 'saddle':
            # half-open circle
        if fp.stability == 'u':
            style = 'wo'
        elif fp.stability == 'c':
            style = 'co'
        else: # 's'
            style = 'ko'

        if flip_coords == True:
            plt.plot(fp.point[y], fp.point[x], style, markersize=markersize, mew=2)
        else:
            plt.plot(fp.point[x], fp.point[y], style, markersize=markersize, mew=2)

def plot_PP_vf_custom(gen, xname, yname, N=20, subdomain=None, scale_exp=0):
    """Draw 2D vector field in (xname, yname) coordinates of given Generator,
    sampling on a uniform grid of n by n points.

    Optional subdomain dictionary specifies axes limits in each variable,
    otherwise Generator's xdomain attribute will be used.

    For systems of dimension > 2, the non-phase plane variables will be held
      constant at their initial condition values set in the Generator.

    Optional scale_exp is an exponent (domain is all reals) which rescales
      size of arrows in case of disparate scales in the vector field. Larger
      values of scale magnify the arrow sizes. For stiff vector fields, values
      from -3 to 3 may be necessary to resolve arrows in certain regions.

    Requires matplotlib 0.99 or later
    """
    assert N > 1
    xdom = gen.xdomain[xname]
    ydom = gen.xdomain[yname]
    if subdomain is not None:
        try:
            xdom = subdomain[xname]
        except KeyError:
            pass
        try:
            ydom = subdomain[yname]
        except KeyError:
            pass
    assert all(dst.isfinite(xdom)), "Must specify a finite domain for x direction"
    assert all(dst.isfinite(ydom)), "Must specify a finite domain for y direction"
    w = xdom[1]-xdom[0]
    h = ydom[1]-ydom[0]

    xdict = gen.initialconditions.copy()

    xix = gen.funcspec.vars.index(xname)
    yix = gen.funcspec.vars.index(yname)

    xs = np.linspace(xdom[0], xdom[1], N)
    ys = np.linspace(ydom[0], ydom[1], N)

    X, Y = np.meshgrid(xs, ys)
    dxs, dys = np.meshgrid(xs, ys)

    dz_big = 0
    vec_dict = {}

    for xi, x in enumerate(xs):
        for yi, y in enumerate(ys):
            xdict.update({xname: x, yname: y})
            dx, dy = gen.Rhs(0, xdict)[[xix, yix]]
            # note order of indices
            dxs[yi,xi] = dx
            dys[yi,xi] = dy
            dz = np.linalg.norm((dx,dy))
            if dz > dz_big:
                dz_big = dz

    plt.quiver(X, Y, dxs, dys, angles='xy', pivot='middle', units='inches',
               scale=dz_big*max(h,w)/(10*np.exp(2*scale_exp)), lw=0.01/np.exp(scale_exp-1),
               headwidth=max(2,1.5/(np.exp(scale_exp-1))),
               #headlength=2*max(2,1.5/(exp(scale_exp-1))),
               width=0.001*max(h,w), minshaft=2, minlength=0.001)

    ax = plt.gca()

    print "xdom: ", xdom
    print "ydom: ", ydom
    ax.set_xlim(xdom)
    ax.set_ylim(ydom)
    plt.draw()

# we must give a name
print "Setting project name..."
DSargs = dst.args(name='M345 - A3 - Bead on a rotating hoop')

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

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

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


# set the domain of integration.
print "Setting domain of integration..."
DSargs.xdomain = {'phi': [-np.pi, np.pi], 'nu': [-4, 4]}
DSargs.tdomain = [0, 50]

DSargs.fnspecs = {'Jacobian': (['t','phi','nu'],
                                """[[0, 1],
                                    [g*cos(phi)*cos(phi) - g*sin(phi)*sin(phi) + cos(phi), -d]]""")}

# an instance of the 'Generator' class.
print "Initializing generator..."
ode_sys  = dst.Generator.Vode_ODEsystem(DSargs)

print "Calculating trajectory..."
traj = ode_sys.compute('Population over time')
print "Extracting data for plotting..."
pts  = traj.sample(dt=0.1)

plt.figure(1)

# phase plane tools are in the Toolbox module


# 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)

# 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-8, max_step=0.1, fps=fp_coords)

# plot the fixed points
for fp in fp_coords:
    fp_ = pp.fixedpoint_2D(ode_sys, dst.Point(fp))
    plot_PP_fps_custom(fp_, 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 this command on your system
plt.show()

Brian, I have fixed up a version of your code that nicely computes the saddle manifolds extending out a long way from the saddle point. It's not fast even using Radau, but I could certainly have better optimized the switch between the two stages of computation with the various parameters after some deeper diagnostics. The ds_perp calculation could also be better automated to use recent curvature information rather than always start from the eigenvector, which will be terribly inaccurate far from the saddle point.

I will post a proper adjustment to the git repo when I have time, but to get you working ASAP I have posted the two files here on dropbox:
https://dl.dropboxusercontent.com/u/2816627/phaseplane.py
https://dl.dropboxusercontent.com/u/2816627/Brian_testing.py
It produces this cool picture: https://dl.dropboxusercontent.com/u/2816627/saddle_manifolds.png
There are several changes to your setup that I've documented in-line using my initials.

I have adjusted the API of the function to make it clearer, and added a lot to the docstring. Please send feedback about it. If I use this in a tutorial, I'll provide a detailed diagram showing how it sets up the problem and what all the parameters control geometrically.

-Rob

Plotting stable and unstable manifolds of a 2D saddles?

Forums

Help

Plotting stable and unstable manifolds of a 2D saddles? document.SUBSCRIPTION_OPTIONS = { "thing": "topic", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Plotting stable and unstable manifolds of a 2D saddles?