In [1]:

```
%matplotlib inline
import random
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt, pi, erf
import seaborn
seaborn.set_context("talk")
import scipy.stats
```

- Know the definition of ordinary differential equation
- Be able to determine the order of an ordinary differential equation
- Convert ordinary differential equations into standard form
- Implement ordinary differential equations into python
- Interpret the output from the ordinary differential equation solver in scipy

The word ordinary means that there are no partial derivative terms nor are there boundary conditions (only initial conditions).

A 1$^\textrm{st}$ order only has first derivatives and a second order has second derivatives and often first derivatives.

*First order Example:*

$\frac{dC}{dt} = -rC$ The rate of chage of concentration of a chemical species as a reaction proceeds

$\frac{dV}{dt} = kt$ Filling a tank at a constant rate ($k$)

$\frac{dV}{dt} = kt - A_{out}\sqrt{2g\frac{V}{A}}$ Filling a tank with a hole in it

*Second order Examples:*

$\frac{d^2x}{dt^2} + k x = 0$ An ideal spring-mass system

$\frac{d^2x}{dt^2} + b \frac{dx}{dt} + k x = 0$ A spring-mass system with dampening

$\frac{d^2T}{dx^2} = x$ Fourier's law of heat conduction

Python can hanlde any order of system of ordinary differential equations. It cannot handle coupled (PDEs) or boundary problems. For those and all other things not in `scipy`

, see this list of numerical methods in python

To solve an ODE in python and all other programming languages/packages:

- Convert to standard form
- Implement the standardized equation as a Python function
- Create a grid of points where you want to evaluate the ODE
- Call
`odeint`

with the function, initial value, and grid

The standard form is:

$$ \frac{dy}{dt} = f(y,t)$$All first order, second order and $n$th order ODEs can be conevrted into this form. Let's see some examples.

The chemical reaction example: $$\frac{dC(t)}{dt} = -rC(t)$$ is already in standard form. So are all the other first-order examples given above.

The trick for the 2nd order form is to create a second dimension, the dimensions of the derivative, which turns the problem into 2D 1st order ODE. Let's see the spring mass system first:

$$\frac{d^2x(t)}{dt^2} + k\, x(t)= 0$$Now we'll use a $x_1(t) = x(t)$ and $x_2(t) = \frac{dx(t)}{dt}$. I'll drop the function notation for simplicity now

$$\frac{dx_2}{dt} + kx_1 = 0$$
$$\frac{dx_2}{dt} = -k x_1$$$$\frac{dx_1}{dt} = x_2 $$

Now we have two dimensions and their ODEs are both in standard form.

Let's see the Fourier equation as another example:

$$\frac{d^2T}{dx^2} = x$$$$T_1 = T$$$$T_2 = \frac{dT}{dx}$$Original diff eq: $$\frac{dT_2}{dx} = x$$

Result of our definition: $$\frac{dT_1}{dx} = T_2$$

We need to implement the right-hand side of this equation:

$$\frac{dC}{dt} = -rC$$In [2]:

```
def rxn_d(C, t, r=1):
#I added a default argument so that I can easily solve
return -r * C
```

Let's see the spring-mass system:

$$\frac{dx_1}{dt} = x_2 $$$$\frac{dx_2}{dt} = -k x_1$$Our function should return the RHS of the governing ODEs for each equation.

In [3]:

```
def ideal_spring(x,t,k=1):
#This will get x as a vector!
return np.array([x[1], x[0] * -k])
```

We'll do Fouier's law now:

$$\frac{dT_1}{dx} = T_2$$$$\frac{dT_2}{dx} = x$$In [4]:

```
def fourier(T, x, k=1):
return np.array([T[1], k * x])
```

Let's solve the reaction system

In [5]:

```
import scipy.integrate
#my initial condition
c_init = 1.0
#where I want to evaluate my function
t_points = np.linspace(0,5, 100)
c_points = scipy.integrate.odeint(rxn_d, c_init, t_points)
plt.plot(t_points, c_points)
plt.xlabel('$t$')
plt.ylabel('$C(t)$')
plt.show()
```

Let's solve the spring-mass system

In [6]:

```
#This is just pasted from above
def ideal_spring(x,t,k=1):
#This will get x as a vector!
return np.array([x[1], x[0] * -k])
x_init = [1.0, 0] #The position is at 1.0 and has no velocity
t_points = np.linspace(0,10, 250)
x_points = scipy.integrate.odeint(ideal_spring, x_init, t_points)
plt.plot(t_points, x_points[:,0]) #Only plot x[0], the actual x-position
plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.show()
```

What about the dampened system?

$$\frac{dx_1}{dt} = x_2 $$$$\frac{dx_2}{dt} = -k x_1 - b x_2$$

In [7]:

```
#This is just pasted from above
def damp_spring(x,t,k=1, b=0.1):
#This will get x as a vector!
return np.array([x[1], x[0] * -k - b * x[1]])
x_init = [1.0, 0] #The position is at 1.0 and has no velocity
t_points = np.linspace(0,100, 250)
x_points = scipy.integrate.odeint(damp_spring, x_init, t_points)
plt.plot(t_points, x_points[:,0]) #Only plot x[0], the actual x-position
plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.show()
```

Since we're in python, we can do some programming to explore the parameter space. Let's do a for-loop to see the effect of different dampenings

In [9]:

```
x_init = [1.0, 0] #The position is at 1.0 and has no velocity
t_points = np.linspace(0,20, 250)
for bi in np.linspace(0.001, 1.0, 5):
x_points = scipy.integrate.odeint(lambda x,t: damp_spring(x,t,k=1.0, b=bi), x_init, t_points)
plt.plot(t_points, x_points[:,0], label='$b_i = {}$'.format(bi)) #Only plot x[0], the actual x-position
plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.ylim(-3, 8) #this is to leave enough room for a legend
plt.legend()
plt.show()
```

We're using the lsoda solver when we call `odeint`

. You can choose between many more by using the more complex `ode`

interface in `scipy`

. We won't cover this in lecture, but this allows you to explore other types of solution methods including the famous runge-kutta method and methods for complex numbers