MAS212 Assignment 2 (A dynamical system): example code


In which we investigate the trajectories of a 2D dynamical system governed by the second-order differential equations $$ \ddot{x} = - \frac{\partial h}{\partial x} , \quad \quad \ddot{y} = - \frac{\partial h}{\partial y} , $$ where $$ h(x,y) = \frac{1}{2} \left(x^2 + y^2 \right) + x^2 y - \frac{1}{3} y^3 . $$

In [1]:
# Part 1 does not require any code and is addressed in the example report.
In [61]:
# Import libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as ode
%matplotlib inline
plt.rcParams.update({'font.size': 14})  # increase the text size.

2(a) A contour plot

In [3]:
## Part 2(a)
xmin, xmax, xpts = -1.5, 1.5, 501
ymin, ymax, ypts = -1.5, 1.5, 501
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
levels = np.linspace(0.0001,1/3,11)   # specify the contour heights using an optional argument of plt.contour()
plt.gcf().set_size_inches(6,6) 
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.title("Contours of $h(x,y)$")
plt.contour(xs.flatten(), ys.flatten(), hpot, levels=levels)
plt.tight_layout()
plt.savefig("hh_contours.pdf")

2(b) First-order form of the equations


Let $p_x = \dot{x}$ and $p_y = \dot{y}$. Taking a derivative with respect to $t$ gives $\dot{p}_x = \ddot{x}$ and $\dot{p}_y = \ddot{y}$. We can replace the left-hand sides of the second-order equations, and rearrange, to get \begin{align} \dot{x} &= p_x , & \dot{p}_x &= -x - 2 x y , \\ \dot{y} &= p_y , & \dot{p}_y &= -y - x^2 + y^2 . \end{align} This is a 4D system of ODEs, with dependent variables $\{x(t), y(t), p_x(t), p_y(t)\}$.

2(c) Example trajectories

In [72]:
from scipy.integrate import odeint

def deq(X, t):
    x, y, px, py = X;  # a 4D system.
    dxdt = px
    dydt = py
    dPxdt = -x - 2*x*y
    dPydt = -y - x**2 + y**2
    return np.array([dxdt, dydt, dPxdt, dPydt])

def energy(X0):
    """Given an initial condition, the function calculates the energy E."""
    x,y,px,py = X0
    E = 1/2*(px**2+py**2) + 1/2*(x**2 + y**2) + x**2*y - 1/3*y**3
    return E

tmax = 100;
tpts = 1000;
ts = np.linspace(0,tmax,tpts)
IC=[0.1,0.0,0.0,0.2]
Xs = odeint(deq, IC, ts)
plt.gcf().set_size_inches((6.5,6))
plt.plot(Xs[:,0], Xs[:,1], "-", linewidth=1.5)
plt.xlabel("x")
plt.ylabel("y")

xmin,xmax=-0.25,0.25
ymin,ymax=-0.25,0.25
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
plt.contour(xs.flatten(), ys.flatten(), hpot, [energy(IC)]);

plt.tight_layout()
plt.savefig("traj1a.pdf")
In [73]:
# Change the initial conditions
IC=[0.1,0,0,0.39991]
Xs = odeint(deq, IC, ts)
plt.gcf().set_size_inches((6.5,6))
plt.plot(Xs[:,0], Xs[:,1], "-", linewidth=1.5)
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.savefig("traj1b.pdf")

xmin,xmax=-0.5,0.5
ymin,ymax=-0.4,0.55
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
plt.contour(xs.flatten(), ys.flatten(), hpot, [energy(IC)]);

plt.tight_layout()
plt.savefig("traj1b.pdf")
In [62]:
tmax = 200;
tpts = 2000;
ts = np.linspace(0,tmax,tpts)
IC=[0.3,0,0,0.45]
Xs = odeint(deq, IC, ts)
plt.gcf().set_size_inches((6.5,6))
plt.plot(Xs[:,0], Xs[:,1], "-", linewidth=1.5)
plt.xlabel("x")
plt.ylabel("y")

xmin,xmax=-0.7,0.7
ymin,ymax=-0.6,0.8
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
plt.contour(xs.flatten(), ys.flatten(), hpot, [energy(IC)]);

plt.tight_layout()
plt.savefig("traj1c.pdf")

3(a) An escape trajectory

In [63]:
# One choice is to use the midpoint method. This is not an essential part of markscheme.

def midpoint(tmax=100, IC=[0.1,0.0,0.0,0.2], h=0.01):
    """My implementation of the midpoint method, described in Lectures 5 and 6."""
    neq = 4
    n = int(tmax/h)
    ts = np.arange(n+1)*h
    Xs = np.zeros((n+1,neq))
    Xs[0,:] = IC
    for k in range(n):
        dXdt = deq(Xs[k,:], ts[k])
        Xhalf = Xs[k,:] + (h/2)*dXdt
        thalf = ts[k] + h/2
        dXdt = deq(Xhalf, thalf)
        Xs[k+1,:] = Xs[k,:] + h*dXdt
    return ts, Xs
In [59]:
IC = [0.2,0.0,0.0,0.56]
print("E/E0 = %.6f" % (6*energy(IC),))
tmax0=39.5
ts, Xs = midpoint(tmax=tmax0, IC=IC, h=0.001)
plt.title("An escape trajectory")
plt.gcf().set_size_inches((6,4.5))
plt.plot(Xs[:,0], Xs[:,1], "-", linewidth=1.5)
plt.xlabel("$x$")
plt.ylabel("$y$")
xmin,xmax=-1.5,1.5
ymin,ymax=-1,1.2
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
plt.contour(xs.flatten(), ys.flatten(), hpot, [energy(IC)]);
plt.tight_layout()
plt.savefig("escape-traj.pdf")
E/E0 = 1.060800
In [64]:
# An alternative choice is to use scipy.integrate.solve_ivp()
# Try using an event condition, as described in the help file for solve_ivp()

from scipy.integrate import solve_ivp

def deq2(t, X):  # The order of arguments is opposite to odeint
    return deq(X,t)

def escaped(t, X): 
    x,y,px,py = X
    return x**2 + y**2 - 1.3**2
escaped.terminal = True

IC = [0.2,0.0,0.0,0.56]
tmax = 100
npts = 1000
ts = np.linspace(0,tmax,npts)
sol = solve_ivp(deq2, (0,tmax), IC, t_eval=ts, events=escaped)

ts = sol.t
Xs = sol.y
plt.gcf().set_size_inches((6.5,6))
plt.plot(Xs[0,:], Xs[1,:], "-", linewidth=1.5)  # Note: columns and rows are switched in the output of solve_ivp() relative to odeint()

plt.title("An escape trajectory")
xmin=-1.0
xmax=1.0
ymin=-0.8
ymax=1.3
xs = np.linspace(xmin,xmax,xpts).reshape((1,xpts))
ys = np.linspace(ymin,ymax,ypts).reshape((ypts,1))
x1s = np.ones(xpts).reshape((1,xpts))
hpot = 1/2*(xs**2 + ys**2) + (xs**2)*ys - 1/3*ys**3*x1s
plt.contour(xs.flatten(), ys.flatten(), hpot, [energy(IC)]);
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.tight_layout()
plt.savefig("escape-traj.pdf")

Escape basins

In [71]:
def escaped(t, X): 
    x,y,px,py = X
    return x**2 + y**2 - 1.3**2
escaped.terminal = True

def testescape(tmax=200.0, IC=[0.1,0.0,0.0,0.2]):
    """This function returns 0,1,2,3, with 0 if the trajectory does not escape,
    and 1 if it escapes 'North', 2 if South-west, 3 if South-east. """
    sol = solve_ivp(deq2, (0,tmax), IC, events=escaped)
    Xs = sol.y
    if np.abs(sol.t[-1] - tmax) < 1e-6:
        trajtype = 0   # still in the basin
    elif Xs[1,-1] >= 1.0:
        trajtype = 1
    elif Xs[0,-1] < 0:
        trajtype = 2
    elif Xs[0,-1] > 0:
        trajtype = 3
    return trajtype
In [57]:
# Warning: this takes a minute or so to calculate
gridpts = 50
xmax, pmax = 0.8, 0.8
x0s = np.linspace(0, xmax, gridpts)
py0s = np.linspace(0, pmax, gridpts)
escapearr = np.zeros((gridpts, gridpts))
tmax = 100.0
for i in range(gridpts):
    print(i)
    for j in range(gridpts):
        trajtype = testescape(tmax=tmax, IC=[x0s[i],0.0,0.0,py0s[j]])
        escapearr[i,j] = trajtype
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
In [70]:
plt.xlabel('$x(0)$')
plt.ylabel('$p_y(0)$')
plt.gcf().set_size_inches(5,6)
plt.imshow(escapearr, cmap='jet', extent=[0, xmax, pmax, 0])
plt.tight_layout()
plt.savefig("escapearr.pdf")