# Worked Solutions for Lab Class 4: Solving ODEs numerically

In [45]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
plt.rcParams.update({'font.size': 14})  # Make the labels larger


### Question 1: First-order Ordinary Differential Equations (ODEs).

In [46]:
# Question 1(a)
def dxdt(x,t):
return -x

h = 0.1;   # h is the spacing in time.
tmax = 5.0; tmin = 0.0;
n = int((tmax-tmin)/h)  # n+1 is the number of points
ts = np.linspace(tmin, tmax, n+1)  # array of t values
x0 = 1.0;    # the initial condition
xexact = np.exp(-ts)  # This is the analytic solution of the ODE
xs = odeint(dxdt, x0, ts)  # The numerical solution.
plt.plot(ts, xs, 'b+');
plt.plot(ts, xexact, 'r-');
plt.xlabel("t"); plt.ylabel("x(t)");

In [47]:
# Now plot the error: the difference between the numerical solution and the analytic solution
# (Here it is necessary to convert 'xs' from a 2D array to a 1D array.)
plt.plot(ts, xs.reshape(n+1) - xexact, 'r-');

In [48]:
# The error in the first plot is of order 10^(-8).
# We can reduce the error by adjusting the optional parameters 'atol' and 'rtol' of odeint.
# These parameters fix the absolute and relative errors that odeint will tolerate.
tol = 1e-11;
xs = odeint(dxdt, x0, ts, atol=tol, rtol=tol);
# Note that the error is now approximately 1000 times smaller than before.
plt.plot(ts, xs.reshape(n+1) - xexact, 'r-');

In [49]:
# Question 1(b).
def dxdt(x,t):
return x*(1-x)
x0s = np.linspace(0, 2, 10)
for x0 in x0s:
xs = odeint(dxdt, x0, ts)
plt.plot(ts, xs)
# x = 1 is a stable equilibrium and x = 0 is the unstable equilibrium
plt.xlabel('t'); plt.ylabel('x(t)');

In [50]:
# Question 1(c).
def dxdt(x,t):
return -1/2*(x-1)*x*(x+1)
x0s = np.linspace(-3, 3, 25)
for x0 in x0s:
xs = odeint(dxdt, x0, ts)
plt.plot(ts, xs)
# x = +1 and x = -1 are stable equilibria, and x = 0 is an unstable equilibrium.
plt.xlabel('t'); plt.ylabel('x(t)');


### Question 2: The Predator-Prey Equations.

In [51]:
# Question 2(a).
a,b,c,d = 1,1,1,1

def predprey(X, t):
# X is an array with two elements. The function should return an array of the derivatives of these elements.
x = X[0]; # rabbits
y = X[1]; # foxes
dxdt =  a*x - b*x*y
dydt = -c*y + d*x*y
return np.array([dxdt, dydt])

h = 0.1; tmax = 20.0; tmin = 0.0;
n = int((tmax-tmin)/h); ts = np.linspace(tmin, tmax, n+1);

ic = (1.0, 1.5)  # initial conditions: (rabbits, foxes).
Xs = odeint(predprey, ic, ts)
# Note that odeint returns a 2D array with 2 columns.
# Rabbits in the first column, foxes in the second.
plt.plot(ts, Xs[:,0], 'r-', label="Rabbits");
plt.plot(ts, Xs[:,1], 'b-', label="Foxes");
plt.xlabel("time"); plt.ylabel("Population"); plt.legend();

In [53]:
# Question 1(b): a phase portrait
r0 = 1.0;  # Initial number of rabbits
nseries = 6;  # number of curves on the phase plot
fs = np.linspace(0.2, 1, nseries)  # Initial number of foxes
for k in range(nseries):
Xs = odeint(predprey, (r0, fs[k]), ts)
plt.plot(Xs[:,0], Xs[:,1]);
plt.xlabel("Rabbits"); plt.ylabel("Foxes");

In [65]:
# Question 1(c): Try varying parameter 'b', the appetite of the foxes
def predprey2(X, t, b=1.0):  # Make b an optional parameter
x = X[0]; y = X[1];
dxdt =  a*x - b*x*y
dydt = -c*y + d*x*y
return np.array([dxdt, dydt])

bs = np.linspace(0.5, 2.0, 4)
ic = (0.5, 1.0)  # Try starting with this initial condition
for bval in bs:
# Pass the value of the 'b' parameter through to predprey2 using the optional argument 'args'
Xs = odeint(predprey2, ic, ts, args=(bval,))
plt.plot(Xs[:,0], Xs[:,1]);
# Changing 'b' changes the position of the stable equilibrium.


### Q3. Second-order equations

Q3(a). After introducing the new variable $y = \dot{x}$, I obtain a 2D set of first-order ODEs: \begin{eqnarray} \dot{x} &=& y \ \dot{y} &=& -2\gamma y - x , \quad x(0) = 1, y(0) = 0. \end{eqnarray}

In [77]:
# Question 3(b)
def oscillator(X, t, gamma = 1.0):
x = X[0]; y = X[1];
dxdt = y
dydt = -2*gamma*y - x
return np.array([dxdt, dydt])

h = 0.1; tmax = 20.0; tmin = 0.0;
n = int((tmax-tmin)/h); ts = np.linspace(tmin, tmax, n+1);
X0 = (1.0, 0.0) # initial conditions
nseries = 10
gammas = np.linspace(0, 2.0, nseries+1)
for g in gammas:
# Pass the parameter 'gamma' via the optional argument 'args'
Xs = odeint(oscillator, X0, ts, args=(g,))
plt.plot(Xs[:,0], Xs[:,1])
plt.axis('equal');
# The point (0,0) is an attractive limit point if gamma > 0.
plt.xlabel('x(t)'); plt.ylabel('y(t)');


### Q4. Van der Pol equation

In [87]:
def vanderPol(X, t, a = 1.0):
x = X[0]; y = X[1];
dxdt = y
dydt = a*(1-x**2)*y - x
return np.array([dxdt, dydt])

h = 0.05; tmax = 50.0; tmin = 0.0;
n = int((tmax-tmin)/h); ts = np.linspace(tmin, tmax, n+1);
X0 = (0.5, 1.0) # initial conditions
nseries = 6
avals = np.linspace(0.2, 2.0, nseries+1)
plt.xlim((-2,2));
for a in avals:
# Pass the parameter 'a' via the optional argument 'args'
Xs = odeint(vanderPol, X0, ts, args=(a,))
plt.plot(Xs[n//2:,0], Xs[n//2:,1])
plt.axis('equal');
# The point (0,0) is an attractive limit point if gamma > 0.
plt.xlabel('x(t)'); plt.ylabel('y(t)');
plt.title("Limit cycles for a = 0.2 to 2.0");


Q4(c). The forced van der Pol equation was investigated in Assignment 2 of 2015/16. See here for a report:
http://sam-dolan.staff.shef.ac.uk/2015/mas212/docs/assign2_report.pdf

In [ ]: