In this notebook we will use Python to solve differential equations numerically.

In [1]:

```
# Our first step is to import any required modules.
# I will use a shortcut to import commonly-used modules in scientific computing:
%pylab inline
```

Let's try a first-order ordinary differential equation (ODE), say: \[\quad \frac{dy}{dx} + y = x, \quad \quad y(0) = 1. \] This has a closed-form solution \[\quad y = x - 1 + 2e^{-x} \] (Exercise: Show this, by first finding the integrating factor.)

We are going to solve this numerically.

First, let's import the "scipy" module and look at the help file for the relevant function, "integrate.odeint",

In [2]:

```
from scipy.integrate import odeint
#odeint? # Uncomment to view the help file for this function
# Define a function which calculates the derivative
def dy_dx(y, x):
return x - y
xs = np.linspace(0,5,100)
y0 = 1.0
ys = odeint(dy_dx, y0, xs)
ys = np.array(ys).flatten()
```

In [3]:

```
plt.plot(xs, ys)
plt.xlabel("x")
plt.ylabel("y");
```

Compare the numerical solution with the analytical solution by showing both on the same plot

In [4]:

```
y_exact = xs - 1 + 2*np.exp(-xs)
y_difference = ys - y_exact
plt.plot(xs, ys, xs, y_exact, "+");
```

Now take a look at the difference between the two series:

In [5]:

```
y_diff = np.abs(y_exact - ys)
plt.semilogy(xs, y_diff)
plt.ylabel("Error")
plt.xlabel("x")
plt.title("Error in numerical integration");
```

Exercise: Experiment with the options of "odeint" to improve the accuracy of the integration.

In [6]:

```
def dU_dx(U, x):
# Here U is a vector such that y=U[0] and z=U[1]. This function should return [y', z']
return [U[1], -2*U[1] - 2*U[0] + np.cos(2*x)]
U0 = [0, 0]
xs = np.linspace(0, 10, 200)
Us = odeint(dU_dx, U0, xs)
ys = Us[:,0]
```

In [7]:

```
plt.plot(xs,ys)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Damped harmonic oscillator");
```

Also known as Lotka-Volterra equations, the predator-prey equations are a pair of first-order non-linear ordinary differential equations. They represent a simplified model of the change in populations of two species which interact via predation. For example, foxes (predators) and rabbits (prey). Let \(x\) and \(y\) represent rabbit and fox populations, respectively. Then

Here \(a\), \(b\), \(c\) and \(d\) are parameters, which are assumed to be positive.

In [8]:

```
a,b,c,d = 1,1,1,1
def dP_dt(P, t):
return [P[0]*(a - b*P[1]), -P[1]*(c - d*P[0])]
ts = np.linspace(0, 12, 100)
P0 = [1.5, 1.0]
Ps = odeint(dP_dt, P0, ts)
prey = Ps[:,0]
predators = Ps[:,1]
```

In [9]:

```
plt.plot(ts, prey, "+", label="Rabbits")
plt.plot(ts, predators, "x", label="Foxes")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend();
```

In [10]:

```
plt.plot(prey, predators, "-")
plt.xlabel("Rabbits")
plt.ylabel("Foxes")
plt.title("Rabbits vs Foxes");
```

In [11]:

```
ic = np.linspace(1.0, 3.0, 21)
for r in ic:
P0 = [r, 1.0]
Ps = odeint(dP_dt, P0, ts)
plt.plot(Ps[:,0], Ps[:,1], "-")
plt.xlabel("Rabbits")
plt.ylabel("Foxes")
plt.title("Rabbits vs Foxes");
```