# phaseplot.py - generate 2D phase portraits

#

# Initial author: Richard M. Murray

# Creation date: 24 July 2011, converted from MATLAB version (2002);

# based on an original version by Kristi Morgansen

"""Generate 2D phase portraits.

This module contains functions for generating 2D phase plots. The base

function for creating phase plane portraits is `~control.phase_plane_plot`,

which generates a phase plane portrait for a 2 state I/O system (with no

inputs). Utility functions are available to customize the individual

elements of a phase plane portrait.

The docstring examples assume the following import commands::

>>> import numpy as np

>>> import control as ct

>>> import control.phaseplot as pp

"""

import math

import warnings

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

from scipy.integrate import odeint

from . import config

from .ctrlplot import ControlPlot, _add_arrows_to_line2D, _get_color, \

_process_ax_keyword, _update_plot_title

from .exception import ControlArgument

from .nlsys import NonlinearIOSystem, find_operating_point, \

input_output_response

__all__ = ['phase_plane_plot', 'phase_plot', 'box_grid']

# Default values for module parameter variables

_phaseplot_defaults = {

'phaseplot.arrows': 2, # number of arrows around curve

'phaseplot.arrow_size': 8, # pixel size for arrows

'phaseplot.arrow_style': None, # set arrow style

'phaseplot.separatrices_radius': 0.1 # initial radius for separatrices

}

def phase_plane_plot(

sys, pointdata=None, timedata=None, gridtype=None, gridspec=None,

plot_streamlines=None, plot_vectorfield=None, plot_streamplot=None,

plot_equilpoints=True, plot_separatrices=True, ax=None,

suppress_warnings=False, title=None, **kwargs

):

"""Plot phase plane diagram.

This function plots phase plane data, including vector fields, stream

lines, equilibrium points, and contour curves.

If none of plot_streamlines, plot_vectorfield, or plot_streamplot are

set, then plot_streamplot is used by default.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot or an array of shape (N, 2)

giving points of at which to plot the vector field.

timedata : int or list of int

Time to simulate each streamline. If a list is given, a different

time can be used for each initial condition in `pointdata`.

gridtype : str, optional

The type of grid to use for generating initial conditions:

'meshgrid' (default) generates a mesh of initial conditions within

the specified boundaries, 'boxgrid' generates initial conditions

along the edges of the boundary, 'circlegrid' generates a circle of

initial conditions around each point in point data.

gridspec : list, optional

If the gridtype is 'meshgrid' and 'boxgrid', `gridspec` gives the

size of the grid in the x and y axes on which to generate points.

If gridtype is 'circlegrid', then `gridspec` is a 2-tuple

specifying the radius and number of points around each point in the

`pointdata` array.

params : dict, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : matplotlib color spec, optional

Plot all elements in the given color (use ``plot_<element>`` =

{'color': c} to set the color in one element of the phase

plot (equilpoints, separatrices, streamlines, etc).

ax : `matplotlib.axes.Axes`, optional

The matplotlib axes to draw the figure on. If not specified and

the current figure has a single axes, that axes is used.

Otherwise, a new figure is created.

Returns

-------

cplt : `ControlPlot` object

Object containing the data that were plotted. See `ControlPlot`

for more detailed information.

cplt.lines : array of list of `matplotlib.lines.Line2D`

Array of list of `matplotlib.artist.Artist` objects:

- lines[0] = list of Line2D objects (streamlines, separatrices).

- lines[1] = Quiver object (vector field arrows).

- lines[2] = list of Line2D objects (equilibrium points).

- lines[3] = StreamplotSet object (lines with arrows).

cplt.axes : 2D array of `matplotlib.axes.Axes`

Axes for each subplot.

cplt.figure : `matplotlib.figure.Figure`

Figure containing the plot.

Other Parameters

----------------

arrows : int

Set the number of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrows']`.

arrow_size : float

Set the size of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrow_size']`.

arrow_style : matplotlib patch

Set the style of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrow_style']`.

dir : str, optional

Direction to draw streamlines: 'forward' to flow forward in time

from the reference points, 'reverse' to flow backward in time, or

'both' to flow both forward and backward. The amount of time to

simulate in each direction is given by the `timedata` argument.

plot_streamlines : bool or dict, optional

If True then plot streamlines based on the pointdata and gridtype.

If set to a dict, pass on the key-value pairs in the dict as

keywords to `streamlines`.

plot_vectorfield : bool or dict, optional

If True then plot the vector field based on the pointdata and

gridtype. If set to a dict, pass on the key-value pairs in the

dict as keywords to `phaseplot.vectorfield`.

plot_streamplot : bool or dict, optional

If True then use `matplotlib.axes.Axes.streamplot` function

to plot the streamlines. If set to a dict, pass on the key-value

pairs in the dict as keywords to `phaseplot.streamplot`.

plot_equilpoints : bool or dict, optional

If True (default) then plot equilibrium points based in the phase

plot boundary. If set to a dict, pass on the key-value pairs in the

dict as keywords to `phaseplot.equilpoints`.

plot_separatrices : bool or dict, optional

If True (default) then plot separatrices starting from each

equilibrium point. If set to a dict, pass on the key-value pairs

in the dict as keywords to `phaseplot.separatrices`.

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.defaults['ctrlplot.rcParams']`.

suppress_warnings : bool, optional

If set to True, suppress warning messages in generating trajectories.

title : str, optional

Set the title of the plot. Defaults to plot type and system name(s).

Notes

-----

The default method for producing streamlines is determined based on which

keywords are specified, with `plot_streamplot` serving as the generic

default. If any of the `arrows`, `arrow_size`, `arrow_style`, or `dir`

keywords are used and neither `plot_streamlines` nor `plot_streamplot` is

set, then `plot_streamlines` will be set to True. If neither

`plot_streamlines` nor `plot_vectorfield` set set to True, then

`plot_streamplot` will be set to True.

"""

# Check for legacy usage of plot_streamlines

streamline_keywords = [

'arrows', 'arrow_size', 'arrow_style', 'dir']

if plot_streamlines is None:

if any([kw in kwargs for kw in streamline_keywords]):

warnings.warn(

"detected streamline keywords; use plot_streamlines to set",

FutureWarning)

plot_streamlines = True

if gridtype not in [None, 'meshgrid']:

warnings.warn(

"streamplots only support gridtype='meshgrid'; "

"falling back to streamlines")

plot_streamlines = True

if plot_streamlines is None and plot_vectorfield is None \

and plot_streamplot is None:

plot_streamplot = True

if plot_streamplot and not plot_streamlines and not plot_vectorfield:

gridspec = gridspec or [25, 25]

# Process arguments

params = kwargs.get('params', None)

sys = _create_system(sys, params)

pointdata = [-1, 1, -1, 1] if pointdata is None else pointdata

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Create axis if needed

user_ax = ax

fig, ax = _process_ax_keyword(user_ax, squeeze=True, rcParams=rcParams)

# Create copy of kwargs for later checking to find unused arguments

initial_kwargs = dict(kwargs)

# Utility function to create keyword arguments

def _create_kwargs(global_kwargs, local_kwargs, **other_kwargs):

new_kwargs = dict(global_kwargs)

new_kwargs.update(other_kwargs)

if isinstance(local_kwargs, dict):

new_kwargs.update(local_kwargs)

return new_kwargs

# Create list for storing outputs

out = np.array([[], None, None, None], dtype=object)

# the maximum zorder of stramlines, vectorfield or streamplot

flow_zorder = None

# Plot out the main elements

if plot_streamlines:

kwargs_local = _create_kwargs(

kwargs, plot_streamlines, gridspec=gridspec, gridtype=gridtype,

ax=ax)

out[0] += streamlines(

sys, pointdata, timedata, _check_kwargs=False,

suppress_warnings=suppress_warnings, **kwargs_local)

new_zorder = max(elem.get_zorder() for elem in out[0])

flow_zorder = max(flow_zorder, new_zorder) if flow_zorder \

else new_zorder

# Get rid of keyword arguments handled by streamlines

for kw in ['arrows', 'arrow_size', 'arrow_style', 'color',

'dir', 'params']:

initial_kwargs.pop(kw, None)

# Reset the gridspec for the remaining commands, if needed

if gridtype not in [None, 'boxgrid', 'meshgrid']:

gridspec = None

if plot_vectorfield:

kwargs_local = _create_kwargs(

kwargs, plot_vectorfield, gridspec=gridspec, ax=ax)

out[1] = vectorfield(

sys, pointdata, _check_kwargs=False, **kwargs_local)

new_zorder = out[1].get_zorder()

flow_zorder = max(flow_zorder, new_zorder) if flow_zorder \

else new_zorder

# Get rid of keyword arguments handled by vectorfield

for kw in ['color', 'params']:

initial_kwargs.pop(kw, None)

if plot_streamplot:

if gridtype not in [None, 'meshgrid']:

raise ValueError(

"gridtype must be 'meshgrid' when using streamplot")

kwargs_local = _create_kwargs(

kwargs, plot_streamplot, gridspec=gridspec, ax=ax)

out[3] = streamplot(

sys, pointdata, _check_kwargs=False, **kwargs_local)

new_zorder = max(out[3].lines.get_zorder(), out[3].arrows.get_zorder())

flow_zorder = max(flow_zorder, new_zorder) if flow_zorder \

else new_zorder

# Get rid of keyword arguments handled by streamplot

for kw in ['color', 'params']:

initial_kwargs.pop(kw, None)

sep_zorder = flow_zorder + 1 if flow_zorder else None

if plot_separatrices:

kwargs_local = _create_kwargs(

kwargs, plot_separatrices, gridspec=gridspec, ax=ax)

kwargs_local['zorder'] = kwargs_local.get('zorder', sep_zorder)

out[0] += separatrices(

sys, pointdata, _check_kwargs=False, **kwargs_local)

sep_zorder = max(elem.get_zorder() for elem in out[0]) if out[0] \

else None

# Get rid of keyword arguments handled by separatrices

for kw in ['arrows', 'arrow_size', 'arrow_style', 'params']:

initial_kwargs.pop(kw, None)

equil_zorder = sep_zorder + 1 if sep_zorder else None

if plot_equilpoints:

kwargs_local = _create_kwargs(

kwargs, plot_equilpoints, gridspec=gridspec, ax=ax)

kwargs_local['zorder'] = kwargs_local.get('zorder', equil_zorder)

out[2] = equilpoints(

sys, pointdata, _check_kwargs=False, **kwargs_local)

# Get rid of keyword arguments handled by equilpoints

for kw in ['params']:

initial_kwargs.pop(kw, None)

# Make sure all keyword arguments were used

if initial_kwargs:

raise TypeError("unrecognized keywords: ", str(initial_kwargs))

if user_ax is None:

if title is None:

title = f"Phase portrait for {sys.name}"

_update_plot_title(title, use_existing=False, rcParams=rcParams)

ax.set_xlabel(sys.state_labels[0])

ax.set_ylabel(sys.state_labels[1])

plt.tight_layout()

return ControlPlot(out, ax, fig)

def vectorfield(

sys, pointdata, gridspec=None, zorder=None, ax=None,

suppress_warnings=False, _check_kwargs=True, **kwargs):

"""Plot a vector field in the phase plane.

This function plots a vector field for a two-dimensional state

space system.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot or an array of shape (N, 2)

giving points of at which to plot the vector field.

gridtype : str, optional

The type of grid to use for generating initial conditions:

'meshgrid' (default) generates a mesh of initial conditions within

the specified boundaries, 'boxgrid' generates initial conditions

along the edges of the boundary, 'circlegrid' generates a circle of

initial conditions around each point in point data.

gridspec : list, optional

If the gridtype is 'meshgrid' and 'boxgrid', `gridspec` gives the

size of the grid in the x and y axes on which to generate points.

If gridtype is 'circlegrid', then `gridspec` is a 2-tuple

specifying the radius and number of points around each point in the

`pointdata` array.

params : dict or list, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : matplotlib color spec, optional

Plot the vector field in the given color.

ax : `matplotlib.axes.Axes`, optional

Use the given axes for the plot, otherwise use the current axes.

Returns

-------

out : Quiver

Other Parameters

----------------

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.defaults['ctrlplot.rcParams']`.

suppress_warnings : bool, optional

If set to True, suppress warning messages in generating trajectories.

zorder : float, optional

Set the zorder for the vectorfield. In not specified, it will be

automatically chosen by `matplotlib.axes.Axes.quiver`.

"""

# Process keywords

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Get system parameters

params = kwargs.pop('params', None)

# Create system from callable, if needed

sys = _create_system(sys, params)

# Determine the points on which to generate the vector field

points, _ = _make_points(pointdata, gridspec, 'meshgrid')

# Create axis if needed

if ax is None:

ax = plt.gca()

# Set the plotting limits

xlim, ylim, maxlim = _set_axis_limits(ax, pointdata)

# Figure out the color to use

color = _get_color(kwargs, ax=ax)

# Make sure all keyword arguments were processed

if _check_kwargs and kwargs:

raise TypeError("unrecognized keywords: ", str(kwargs))

# Generate phase plane (quiver) data

vfdata = np.zeros((points.shape[0], 4))

sys._update_params(params)

for i, x in enumerate(points):

vfdata[i, :2] = x

vfdata[i, 2:] = sys._rhs(0, x, np.zeros(sys.ninputs))

with plt.rc_context(rcParams):

out = ax.quiver(

vfdata[:, 0], vfdata[:, 1], vfdata[:, 2], vfdata[:, 3],

angles='xy', color=color, zorder=zorder)

return out

def streamplot(

sys, pointdata, gridspec=None, zorder=None, ax=None, vary_color=False,

vary_linewidth=False, cmap=None, norm=None, suppress_warnings=False,

_check_kwargs=True, **kwargs):

"""Plot streamlines in the phase plane.

This function plots the streamlines for a two-dimensional state

space system using the `matplotlib.axes.Axes.streamplot` function.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot.

gridspec : list, optional

Specifies the size of the grid in the x and y axes on which to

generate points.

params : dict or list, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : matplotlib color spec, optional

Plot the vector field in the given color.

ax : `matplotlib.axes.Axes`, optional

Use the given axes for the plot, otherwise use the current axes.

Returns

-------

out : StreamplotSet

Containter object with lines and arrows contained in the

streamplot. See `matplotlib.axes.Axes.streamplot` for details.

Other Parameters

----------------

cmap : str or Colormap, optional

Colormap to use for varying the color of the streamlines.

norm : `matplotlib.colors.Normalize`, optional

Normalization map to use for scaling the colormap and linewidths.

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.default['ctrlplot.rcParams']`.

suppress_warnings : bool, optional

If set to True, suppress warning messages in generating trajectories.

vary_color : bool, optional

If set to True, vary the color of the streamlines based on the

magnitude of the vector field.

vary_linewidth : bool, optional.

If set to True, vary the linewidth of the streamlines based on the

magnitude of the vector field.

zorder : float, optional

Set the zorder for the streamlines. In not specified, it will be

automatically chosen by `matplotlib.axes.Axes.streamplot`.

"""

# Process keywords

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Get system parameters

params = kwargs.pop('params', None)

# Create system from callable, if needed

sys = _create_system(sys, params)

# Determine the points on which to generate the streamplot field

points, gridspec = _make_points(pointdata, gridspec, 'meshgrid')

grid_arr_shape = gridspec[::-1]

xs = points[:, 0].reshape(grid_arr_shape)

ys = points[:, 1].reshape(grid_arr_shape)

# Create axis if needed

if ax is None:

ax = plt.gca()

# Set the plotting limits

xlim, ylim, maxlim = _set_axis_limits(ax, pointdata)

# Figure out the color to use

color = _get_color(kwargs, ax=ax)

# Make sure all keyword arguments were processed

if _check_kwargs and kwargs:

raise TypeError("unrecognized keywords: ", str(kwargs))

# Generate phase plane (quiver) data

sys._update_params(params)

us_flat, vs_flat = np.transpose(

[sys._rhs(0, x, np.zeros(sys.ninputs)) for x in points])

us, vs = us_flat.reshape(grid_arr_shape), vs_flat.reshape(grid_arr_shape)

magnitudes = np.linalg.norm([us, vs], axis=0)

norm = norm or mpl.colors.Normalize()

normalized = norm(magnitudes)

cmap = plt.get_cmap(cmap)

with plt.rc_context(rcParams):

default_lw = plt.rcParams['lines.linewidth']

min_lw, max_lw = 0.25*default_lw, 2*default_lw

linewidths = normalized * (max_lw - min_lw) + min_lw \

if vary_linewidth else None

color = magnitudes if vary_color else color

out = ax.streamplot(

xs, ys, us, vs, color=color, linewidth=linewidths, cmap=cmap,

norm=norm, zorder=zorder)

return out

def streamlines(

sys, pointdata, timedata=1, gridspec=None, gridtype=None, dir=None,

zorder=None, ax=None, _check_kwargs=True, suppress_warnings=False,

**kwargs):

"""Plot stream lines in the phase plane.

This function plots stream lines for a two-dimensional state space

system.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot or an array of shape (N, 2)

giving points of at which to plot the vector field.

timedata : int or list of int

Time to simulate each streamline. If a list is given, a different

time can be used for each initial condition in `pointdata`.

gridtype : str, optional

The type of grid to use for generating initial conditions:

'meshgrid' (default) generates a mesh of initial conditions within

the specified boundaries, 'boxgrid' generates initial conditions

along the edges of the boundary, 'circlegrid' generates a circle of

initial conditions around each point in point data.

gridspec : list, optional

If the gridtype is 'meshgrid' and 'boxgrid', `gridspec` gives the

size of the grid in the x and y axes on which to generate points.

If gridtype is 'circlegrid', then `gridspec` is a 2-tuple

specifying the radius and number of points around each point in the

`pointdata` array.

dir : str, optional

Direction to draw streamlines: 'forward' to flow forward in time

from the reference points, 'reverse' to flow backward in time, or

'both' to flow both forward and backward. The amount of time to

simulate in each direction is given by the `timedata` argument.

params : dict or list, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : str

Plot the streamlines in the given color.

ax : `matplotlib.axes.Axes`, optional

Use the given axes for the plot, otherwise use the current axes.

Returns

-------

out : list of Line2D objects

Other Parameters

----------------

arrows : int

Set the number of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrows']`.

arrow_size : float

Set the size of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrow_size']`.

arrow_style : matplotlib patch

Set the style of arrows to plot along the streamlines. The default

value can be set in `config.defaults['phaseplot.arrow_style']`.

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.defaults['ctrlplot.rcParams']`.

suppress_warnings : bool, optional

If set to True, suppress warning messages in generating trajectories.

zorder : float, optional

Set the zorder for the streamlines. In not specified, it will be

automatically chosen by `matplotlib.axes.Axes.plot`.

"""

# Process keywords

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Get system parameters

params = kwargs.pop('params', None)

# Create system from callable, if needed

sys = _create_system(sys, params)

# Parse the arrows keyword

arrow_pos, arrow_style = _parse_arrow_keywords(kwargs)

# Determine the points on which to generate the streamlines

points, gridspec = _make_points(pointdata, gridspec, gridtype=gridtype)

if dir is None:

dir = 'both' if gridtype == 'meshgrid' else 'forward'

# Create axis if needed

if ax is None:

ax = plt.gca()

# Set the axis limits

xlim, ylim, maxlim = _set_axis_limits(ax, pointdata)

# Figure out the color to use

color = _get_color(kwargs, ax=ax)

# Make sure all keyword arguments were processed

if _check_kwargs and kwargs:

raise TypeError("unrecognized keywords: ", str(kwargs))

# Create reverse time system, if needed

if dir != 'forward':

revsys = NonlinearIOSystem(

lambda t, x, u, params: -np.asarray(sys.updfcn(t, x, u, params)),

sys.outfcn, states=sys.nstates, inputs=sys.ninputs,

outputs=sys.noutputs, params=sys.params)

else:

revsys = None

# Generate phase plane (streamline) data

out = []

for i, X0 in enumerate(points):

# Create the trajectory for this point

timepts = _make_timepts(timedata, i)

traj = _create_trajectory(

sys, revsys, timepts, X0, params, dir,

gridtype=gridtype, gridspec=gridspec, xlim=xlim, ylim=ylim,

suppress_warnings=suppress_warnings)

# Plot the trajectory (if there is one)

if traj.shape[1] > 1:

with plt.rc_context(rcParams):

out += ax.plot(traj[0], traj[1], color=color, zorder=zorder)

# Add arrows to the lines at specified intervals

_add_arrows_to_line2D(

ax, out[-1], arrow_pos, arrowstyle=arrow_style, dir=1)

return out

def equilpoints(

sys, pointdata, gridspec=None, color='k', zorder=None, ax=None,

_check_kwargs=True, **kwargs):

"""Plot equilibrium points in the phase plane.

This function plots the equilibrium points for a planar dynamical system.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot or an array of shape (N, 2)

giving points of at which to plot the vector field.

gridtype : str, optional

The type of grid to use for generating initial conditions:

'meshgrid' (default) generates a mesh of initial conditions within

the specified boundaries, 'boxgrid' generates initial conditions

along the edges of the boundary, 'circlegrid' generates a circle of

initial conditions around each point in point data.

gridspec : list, optional

If the gridtype is 'meshgrid' and 'boxgrid', `gridspec` gives the

size of the grid in the x and y axes on which to generate points.

If gridtype is 'circlegrid', then `gridspec` is a 2-tuple

specifying the radius and number of points around each point in the

`pointdata` array.

params : dict or list, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : str

Plot the equilibrium points in the given color.

ax : `matplotlib.axes.Axes`, optional

Use the given axes for the plot, otherwise use the current axes.

Returns

-------

out : list of Line2D objects

Other Parameters

----------------

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.defaults['ctrlplot.rcParams']`.

zorder : float, optional

Set the zorder for the equilibrium points. In not specified, it will

be automatically chosen by `matplotlib.axes.Axes.plot`.

"""

# Process keywords

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Get system parameters

params = kwargs.pop('params', None)

# Create system from callable, if needed

sys = _create_system(sys, params)

# Create axis if needed

if ax is None:

ax = plt.gca()

# Set the axis limits

xlim, ylim, maxlim = _set_axis_limits(ax, pointdata)

# Determine the points on which to generate the vector field

gridspec = [5, 5] if gridspec is None else gridspec

points, _ = _make_points(pointdata, gridspec, 'meshgrid')

# Make sure all keyword arguments were processed

if _check_kwargs and kwargs:

raise TypeError("unrecognized keywords: ", str(kwargs))

# Search for equilibrium points

equilpts = _find_equilpts(sys, points, params=params)

# Plot the equilibrium points

out = []

for xeq in equilpts:

with plt.rc_context(rcParams):

out += ax.plot(

xeq[0], xeq[1], marker='o', color=color, zorder=zorder)

return out

def separatrices(

sys, pointdata, timedata=None, gridspec=None, zorder=None, ax=None,

_check_kwargs=True, suppress_warnings=False, **kwargs):

"""Plot separatrices in the phase plane.

This function plots separatrices for a two-dimensional state space

system.

Parameters

----------

sys : `NonlinearIOSystem` or callable(t, x, ...)

I/O system or function used to generate phase plane data. If a

function is given, the remaining arguments are drawn from the

`params` keyword.

pointdata : list or 2D array

List of the form [xmin, xmax, ymin, ymax] describing the

boundaries of the phase plot or an array of shape (N, 2)

giving points of at which to plot the vector field.

timedata : int or list of int

Time to simulate each streamline. If a list is given, a different

time can be used for each initial condition in `pointdata`.

gridtype : str, optional

The type of grid to use for generating initial conditions:

'meshgrid' (default) generates a mesh of initial conditions within

the specified boundaries, 'boxgrid' generates initial conditions

along the edges of the boundary, 'circlegrid' generates a circle of

initial conditions around each point in point data.

gridspec : list, optional

If the gridtype is 'meshgrid' and 'boxgrid', `gridspec` gives the

size of the grid in the x and y axes on which to generate points.

If gridtype is 'circlegrid', then `gridspec` is a 2-tuple

specifying the radius and number of points around each point in the

`pointdata` array.

params : dict or list, optional

Parameters to pass to system. For an I/O system, `params` should be

a dict of parameters and values. For a callable, `params` should be

dict with key 'args' and value given by a tuple (passed to callable).

color : matplotlib color spec, optional

Plot the separatrices in the given color. If a single color

specification is given, this is used for both stable and unstable

separatrices. If a tuple is given, the first element is used as

the color specification for stable separatrices and the second

element for unstable separatrices.

ax : `matplotlib.axes.Axes`, optional

Use the given axes for the plot, otherwise use the current axes.

Returns

-------

out : list of Line2D objects

Other Parameters

----------------

rcParams : dict

Override the default parameters used for generating plots.

Default is set by `config.defaults['ctrlplot.rcParams']`.

suppress_warnings : bool, optional

If set to True, suppress warning messages in generating trajectories.

zorder : float, optional

Set the zorder for the separatrices. In not specified, it will be

automatically chosen by `matplotlib.axes.Axes.plot`.

Notes

-----

The value of `config.defaults['separatrices_radius']` is used to set the

offset from the equilibrium point to the starting point of the separatix

traces, in the direction of the eigenvectors evaluated at that

equilibrium point.

"""

# Process keywords

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

# Get system parameters

params = kwargs.pop('params', None)

# Create system from callable, if needed

sys = _create_system(sys, params)

# Parse the arrows keyword

arrow_pos, arrow_style = _parse_arrow_keywords(kwargs)

# Determine the initial states to use in searching for equilibrium points

gridspec = [5, 5] if gridspec is None else gridspec

points, _ = _make_points(pointdata, gridspec, 'meshgrid')

# Find the equilibrium points

equilpts = _find_equilpts(sys, points, params=params)

radius = config._get_param('phaseplot', 'separatrices_radius')

# Create axis if needed

if ax is None:

ax = plt.gca()

# Set the axis limits

xlim, ylim, maxlim = _set_axis_limits(ax, pointdata)

# Figure out the color to use for stable, unstable subspaces

color = _get_color(kwargs)

match color:

case None:

stable_color = 'r'

unstable_color = 'b'

case (stable_color, unstable_color) | [stable_color, unstable_color]:

pass

case single_color:

stable_color = unstable_color = single_color

# Make sure all keyword arguments were processed

if _check_kwargs and kwargs:

raise TypeError("unrecognized keywords: ", str(kwargs))

# Create a "reverse time" system to use for simulation

revsys = NonlinearIOSystem(

lambda t, x, u, params: -np.array(sys.updfcn(t, x, u, params)),

sys.outfcn, states=sys.nstates, inputs=sys.ninputs,

outputs=sys.noutputs, params=sys.params)

# Plot separatrices by flowing backwards in time along eigenspaces

out = []

for i, xeq in enumerate(equilpts):

# Figure out the linearization and eigenvectors

evals, evecs = np.linalg.eig(sys.linearize(xeq, 0, params=params).A)

# See if we have real eigenvalues (=> evecs are meaningful)

if evals[0].imag > 0:

continue

# Create default list of time points

if timedata is not None:

timepts = _make_timepts(timedata, i)

# Generate the traces

for j, dir in enumerate(evecs.T):

# Figure out time vector if not yet computed

if timedata is None:

timescale = math.log(maxlim / radius) / abs(evals[j].real)

timepts = np.linspace(0, timescale)

# Run the trajectory starting in eigenvector directions

for eps in [-radius, radius]:

x0 = xeq + dir * eps

if evals[j].real < 0:

traj = _create_trajectory(

sys, revsys, timepts, x0, params, 'reverse',

gridtype='boxgrid', xlim=xlim, ylim=ylim,

suppress_warnings=suppress_warnings)

color = stable_color

linestyle = '--'

elif evals[j].real > 0:

traj = _create_trajectory(

sys, revsys, timepts, x0, params, 'forward',

gridtype='boxgrid', xlim=xlim, ylim=ylim,

suppress_warnings=suppress_warnings)

color = unstable_color

linestyle = '-'

# Plot the trajectory (if there is one)

if traj.shape[1] > 1:

with plt.rc_context(rcParams):

out += ax.plot(

traj[0], traj[1], color=color,

linestyle=linestyle, zorder=zorder)

# Add arrows to the lines at specified intervals

with plt.rc_context(rcParams):

_add_arrows_to_line2D(

ax, out[-1], arrow_pos, arrowstyle=arrow_style,

dir=1)

return out

#

# User accessible utility functions

#

# Utility function to generate boxgrid (in the form needed here)

def boxgrid(xvals, yvals):

"""Generate list of points along the edge of box.

points = boxgrid(xvals, yvals) generates a list of points that

corresponds to a grid given by the cross product of the x and y values.

Parameters

----------

xvals, yvals : 1D array_like

Array of points defining the points on the lower and left edges of

the box.

Returns

-------

grid : 2D array

Array with shape (p, 2) defining the points along the edges of the

box, where p is the number of points around the edge.

"""

return np.array(

[(x, yvals[0]) for x in xvals[:-1]] + # lower edge

[(xvals[-1], y) for y in yvals[:-1]] + # right edge

[(x, yvals[-1]) for x in xvals[:0:-1]] + # upper edge

[(xvals[0], y) for y in yvals[:0:-1]] # left edge

)

# Utility function to generate meshgrid (in the form needed here)

# TODO: add examples of using grid functions directly

def meshgrid(xvals, yvals):

"""Generate list of points forming a mesh.

points = meshgrid(xvals, yvals) generates a list of points that

corresponds to a grid given by the cross product of the x and y values.

Parameters

----------

xvals, yvals : 1D array_like

Array of points defining the points on the lower and left edges of

the box.

Returns

-------

grid : 2D array

Array of points with shape (n * m, 2) defining the mesh.

"""

xvals, yvals = np.meshgrid(xvals, yvals)

grid = np.zeros((xvals.shape[0] * xvals.shape[1], 2))

grid[:, 0] = xvals.reshape(-1)

grid[:, 1] = yvals.reshape(-1)

return grid

# Utility function to generate circular grid

def circlegrid(centers, radius, num):

"""Generate list of points around a circle.

points = circlegrid(centers, radius, num) generates a list of points

that form a circle around a list of centers.

Parameters

----------

centers : 2D array_like

Array of points with shape (p, 2) defining centers of the circles.

radius : float

Radius of the points to be generated around each center.

num : int

Number of points to generate around the circle.

Returns

-------

grid : 2D array

Array of points with shape (p * num, 2) defining the circles.