Generating a phase portrait of a linear ODE

readme.md

Comments on script structure

1. Setting functions for solving ODE.
2. Solve for few conditions and plot phase trajectories.
3. Visualize the vector field.
4. Visualize eigenvectors (useful for saddles and nodes)

Setting the function for ODE solver

Here we use odeint, which requires function of (X,t) arguments, where X - input vector. We use auxilary functions P and Q to make the code human-readable.

P = lambda x,y: -3*x + 2*y
Q = lambda x,y: +2*x - 4*y
fun = lambda p,t: [ P(p[0],p[1]), Q(p[0],p[1]) ]

Of course, we can write fun directly:

fun = lambda p,t: [ -3 * p[0] + 2 * p[1],
                    +2 * p[0] - 4 * p[1] ]

You can use whichever form you like best.

Solving an equation system & draw the solution

Our solver odeint takes function, inital state vector (Mx1) and time array (Nx1), where M - our phase space dimensionality, N - number of our timepoints. It returns (NxM) matrix, where every column is M-th variable evolution over N timepoints.

In our case M=2.

tt = np.linspace(0,100, 1000)
xy0 = [1,0]
zz = odeint(fun, xy0, tt)

To draw the phase trajectory:

plt.plot(zz[:,0], zz[:,1])
plt.plot([xy0[0]], [xy0[1]], 'o', mec='k', ms=4)

I highly recommend to put some marker into the starting point to make trajectory orientation visible.

To draw the kinetics:

plt.plot(tt, zz[:,0])
plt.plot(tt, zz[:,1])

Vector field drawing

We draw the field using quiver function. It draw arrows for given points and take arguments in the following order:

• array of X coordinates of the points
• array of Y coordinates of the points
• array of arrow vector X coordinates
• array of arrow vector Y coordinates

First, we define a grid with meshgrid, which takes X array (X1, X2, .. Xn) and Y array Y1, Y2, .. Yn) and returns an array of X and Y coordinate of every point in the grid (X1, Y1) (X1, Y2) ... (X1, Yn) (X2, Y1) (X2, Y2) ...

xx = np.linspace(-5,5,20)
yy = np.linspace(-3,4,20)
Xm,Ym = np.meshgrid(xx,yy)

Next, we feed grid coordinates Xm,Ym to quiver as point coordinates and P(Xm,Ym), Q(Xm,Ym) as arrow vector components:

plt.quiver(Xm,Ym, P(Xm,Ym), Q(Xm,Ym), scale_units='xy', angles='xy')

Note that using scale_units='xy', angles='xy' is very important for correct arrow orientation at the phase portrait.

Plotting eigenvectors

Eigenvectors of the Jacoby matrix are very useful for plotting some objects in a phase space, e.g., separatrix manifold.

We find these vectors by executing eig function, which return both eigenvalues and eigenvectors. The latter are in matrix format in columns (!).

Here, we plot 1st vector vecs[:,0] and print the corresponding eigenvalue vals[0]:

vals, vecs = np.linalg.eig(mtx)
plt.plot([ 0, vecs[0,0] ] , [ 0, vecs[1,0] ], 
      '--r', label='$\lambda_1$ = %.3f' % vals[0])

draw_node.py

#!/usr/bin/env python3

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

# Our equation system has a form:
# x' = P(x,y)
# y' = Q(x,y)

P = lambda x,y: -3*x + 2*y
Q = lambda x,y: +2*x - 4*y


fun = lambda p,t: [ P(p[0],p[1]), Q(p[0],p[1]) ]
#            p - vector input [x,y] 
#            : - output vector [x', y']

plt.figure(figsize=(7,7))

# plot few phase trajectories
x0_ = 1.0
tt = np.linspace(0,100, 1000)
for y0_ in [-1,0,1,2,3,4]:
    xy0 = [x0_, y0_]
    zz = odeint(fun, xy0, tt)
    plt.plot(zz[:,0], zz[:,1])
    plt.plot([xy0[0]], [xy0[1]], 'o', mec='k', ms=4)

# visualize the vector field
xx = np.linspace(-5,5,20)
yy = np.linspace(-3,4,20)
Xm,Ym = np.meshgrid(xx,yy)
plt.quiver(Xm,Ym, P(Xm,Ym), Q(Xm,Ym), scale_units='xy', angles='xy')

# visualize eigenvectors of the Jacoby matrix
mtx = np.array([[-3, 2],[2,-4]])
vals, vecs = np.linalg.eig(mtx)
plt.plot([vecs[0,0]*-3, vecs[0,0]*3] , [vecs[1,0]*-3, vecs[1,0]*3], 
             '--r', label='$\lambda_1$ = %.3f' % vals[0])
plt.plot([vecs[0,1]*-3, vecs[0,1]*3] , [vecs[1,1]*-3, vecs[1,1]*3], 
             '--b', label='$\lambda_2$ = %.3f' % vals[1])

# decorations
plt.plot([0,0],[-10,10], 'k-', lw=2)
plt.plot([-10,10],[0,0], 'k-', lw=2)
plt.legend()
plt.xlim(-5,5)
plt.ylim(-3,4)
plt.xlabel('X')
plt.ylabel('Y')

plt.show()

node.png