Phase Portraits¶

In this notebook we’ll look at how to generate phase portraits. A phase diagram shows the trajectories that a dynamical system can take through its phase space. For any system that obeys the markov property we can construct such a diagram, with one dimension for each of the system’s stocks.

Ingredients¶

Libraries¶

In this analysis we’ll use:

numpy to create the grid of points that we’ll sample over
matplotlib via ipython’s pylab magic to construct the quiver plot

%pylab inline
import pysd
import numpy as np

Populating the interactive namespace from numpy and matplotlib

Example Model¶

A single pendulum makes for a great example of oscillatory behavior. The following model constructs equations of motion for the pendulum in radial coordinates, simplifying the model from one constructed in cartesian coordinates. In this model, the two state elements are the angle of the pendulum and the derivative of that angle, its angular velocity.

model = pysd.read_vensim('../../models/Pendulum/Single_Pendulum.mdl')

Recipe¶

Step 1: Define the range over which to plot¶

We’ve defined the angle \(0\) to imply that the pendulum is hanging straight down. To show the full range of motion, we’ll show the angles ragning from \(-1.5\pi\) to \(+1.5\pi\), which will allow us to see the the points of stable and unstable equilibrium.

We also want to sample over a range of angular velocities that are reasonable for our system, for the given values of length and weight.

angular_position = np.linspace(-1.5*np.pi, 1.5*np.pi, 60)
angular_velocity = np.linspace(-2, 2, 20)

Numpy’s meshgrid lets us construct a 2d sample space based upon our arrays

apv, avv = np.meshgrid(angular_position, angular_velocity)

Step 2: Calculate the state space derivatives at a point¶

We’ll define a helper function, which given a point in the state space, will tell us what the derivatives of the state elements will be. One way to do this is to run the model over a single timestep, and extract the derivative information. In this case, the model’s stocks have only one inflow/outflow, so this is the derivative value.

As the derivative of the angular position is just the angular velocity, whatever we pass in for the av parameter should be returned to us as the derivative of ap.

def derivatives(ap, av):
    ret = model.run(params={'angular_position': ap,
                            'angular_velocity': av},
                    return_timestamps=[0,1],
                    return_columns=['change_in_angular_position',
                                    'change_in_angular_velocity'])

    return tuple(ret.loc[0].values)

derivatives(0,1)

(1.0, -0.0)

Step 3: Calculate the state space derivatives across our sample space¶

We can use numpy’s vectorize to make the function accept the 2d sample space we have just created. Now we can generate the derivative of angular position vector dapv and that of the angular velocity vector davv. As before, the derivative of the angular posiiton should be equal to the angular velocity. We check that the vectors are equal.

vderivatives = np.vectorize(derivatives)

dapv, davv = vderivatives(apv, avv)
(dapv == avv).all()

True

Step 4: Plot the phase portrait¶

Now we have everything we need to draw the phase portrait. We’ll use matplotlib’s quiver function, which wants as arguments the grid of x and y coordinates, and the derivatives of these coordinates.

In the plot we see the locations of stable and unstable equilibria, and can eyeball the trajectories that the system will take through the state space by following the arrows.

plt.figure(figsize=(18,6))
plt.quiver(apv, avv, dapv, davv, color='b', alpha=.75)
plt.box('off')
plt.xlim(-1.6*np.pi, 1.6*np.pi)
plt.xlabel('Radians', fontsize=14)
plt.ylabel('Radians/Second', fontsize=14)
plt.title('Phase portrait for a simple pendulum', fontsize=16);