# Lecture 9

## Differential Equations I:

### Separation of Variables


```python
import numpy as np
import sympy as sp
import scipy.integrate
sp.init_printing()
##################################################
##### Matplotlib boilerplate for consistency #####
##################################################
from ipywidgets import interact
from ipywidgets import FloatSlider
from matplotlib import pyplot as plt

%matplotlib inline

from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg')

global_fig_width = 10
global_fig_height = global_fig_width / 1.61803399
font_size = 12

plt.rcParams['axes.axisbelow'] = True
plt.rcParams['axes.edgecolor'] = '0.8'
plt.rcParams['axes.grid'] = True
plt.rcParams['axes.labelpad'] = 8
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['axes.titlepad'] = 16.0
plt.rcParams['axes.titlesize'] = font_size * 1.4
plt.rcParams['figure.figsize'] = (global_fig_width, global_fig_height)
plt.rcParams['font.sans-serif'] = ['Computer Modern Sans Serif', 'DejaVu Sans', 'sans-serif']
plt.rcParams['font.size'] = font_size
plt.rcParams['grid.color'] = '0.8'
plt.rcParams['grid.linestyle'] = 'dashed'
plt.rcParams['grid.linewidth'] = 2
plt.rcParams['lines.dash_capstyle'] = 'round'
plt.rcParams['lines.dashed_pattern'] = [1, 4]
plt.rcParams['xtick.labelsize'] = font_size
plt.rcParams['xtick.major.pad'] = 4
plt.rcParams['xtick.major.size'] = 0
plt.rcParams['ytick.labelsize'] = font_size
plt.rcParams['ytick.major.pad'] = 4
plt.rcParams['ytick.major.size'] = 0
##################################################
```

## What is a Differential Equation?

If $y$ is a variable that depends on the variable $x$, then a **differential equation in $x$ and $y$** is one that relates $\displaystyle \frac{{\rm d}y}{{\rm d}x}$ to $x$, to $y$, or to both.

We write
$$\frac{{\rm d}y}{{\rm d}x}=f(x,y).$$

## How do we Find the Solutions to a Differential Equation?

The simplest differential equations are very similar to integration.

$$\frac{{\rm d}y}{{\rm d}x} = x^3 + 5x^2 + 1$$

Since all the terms in $y$ (just the derivative in this case) on the left, and all the terms involving $x$ on the right.

$${\rm d}y = (x^3 + 5x^2 + 1) ~{\rm d}x$$


Then we integrate both sides
$$\int {\rm d}y = \int (x^3 + 5x^2 + 1) ~{\rm d}x$$


The general solution to this equation is
$$y + A = \int x^3 + 5x^2 + 1\,{\rm d}x = {1\over 4}x^4 + {5\over 3}x^3 + x + B$$
where $A$ and $B$ are a arbitrary constants.

If we are given the **initial conditions**, say, $y=0$ when $x=0$ (often written as $y(0)=0$), 
we can determine the value for the constants; in this example $A=B$.

This approach can be generalised into what is known as the **Separation of Variables Method**.
---------------------------------------------------------------------------------------------

## Example 1:
$${{\rm d}y\over {\rm d}x} = 2x + c$$

The equation is already separated, with all terms in $y$ (just the derivative in 
this case) on the left, and terms involving $x$ on the right.


We then integrate both sides with respect to $x$:

$$ \int \frac{{\rm d}y}{{\rm d}x}\,{\rm d}x = \int (2x + c)\, {\rm d}x $$

$$y + k_1 = x^2 + cx + k_2$$

$$y = x^2 + cx + A\qquad\mbox{with}\qquad A=k_2-k_1$$

### Initial conditions

$$y = x^2 + cx + A$$

If we are told that $y(1) = c$, can we calculate $A$?

We have that

$$y = x^2 + cx + A\;\qquad\mbox{with}\qquad y(1) = c$$

so:

$$c = 1^2 + c\cdot 1 + A$$

$$A = -1$$

Is this right?


```python
c = 1.0
A = -1.0
x0 = 1.0
x1 = 10.0
x = np.linspace(x0,x1,100)
y_exact = x**2 + c*x + A

def dydx(y, x):
    return 2*x + c
y_numerical = scipy.integrate.odeint(dydx, c, x)

```


```python
plt.plot(x,y_exact,lw=2, c='r', label='exact')
plt.plot(x,y_numerical,lw=2, c='k', ls='--',label='numerical')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

```

## Example 2:

$${{\rm d}y\over {\rm d}x} = 5 x^{{3\over 2}} y $$

Separate the variables so that all the terms in $x$ are on one side of the equation and all the terms in $y$ are on the other side:

$${1\over y}\frac{{\rm d}y}{{\rm d}x} = 5 x^{{3\over 2}} $$


$${1\over y}\frac{{\rm d}y}{{\rm d}x} = 5 x^{{3\over 2}} $$

Integrating both sides with respect to $x$:

$$\int {1\over y}\frac{{\rm d}y}{{\rm d}x}\, {\rm d}x = \int 5x^{{3\over 2}}\,{\rm d}x \quad
  \Rightarrow\quad \int {1\over y}\,{\rm d}y = \int 5x^{{3\over 2}}\, {\rm d}x$$

Now integrate on each side:
$${\ln y + A} = { 5x^{{5\over 2}}\over 5/2} + B = 2x^{{5\over 2}} + B$$

So:

$$\ln y = 2x^{{5\over 2}} + C{\rm\quad where\quad}C=B-A$$

## Example 3:

$$e^x \frac{{\rm d}y}{{\rm d}x} = {x \over y}$$


Rearrange to get:

$$\int y\,{\rm d}y = \int x e^{-x}\,{\rm d}x$$

We have

$$\int y\,{\rm d}y = \int x e^{-x}\,{\rm d}x$$

Now integrate by parts, $\displaystyle\biggl(\int uv'\,{\rm d}x = [uv]-\int u'v\,{\rm d}x\biggr)$,

with $u=x$ and $v'=e^{-x}$:

\begin{align*}
{y^2\over 2} + c_1 &= -x e ^{-x} + \int e ^{-x}\,{\rm d}x \\\\
      &= -x e ^{-x} - e ^{-x} + c_2 \\\\
  y^2 &= -2 e ^{-x}(1+x) + c_3 \\\\
  y &= \pm\sqrt{-2 e^{-x} (1+x) + c_3}
\end{align*}

Let's check this answer: $y = \pm\sqrt{-2 e^{-x} (1+x) + c_3}$.

Substitute in initial conditions $y(0) = 0$ gives $c_3 = 2$


```python
c3 = 2.0
x0 = 0.0
x1 = 1.0
x = np.linspace(x0,x1,100)

y_exact = np.sqrt(-2*np.exp(-x)*(1 + x) + c3)

def dydx(y, x):
    return np.exp(-x) * x / y

y_numerical = scipy.integrate.odeint(dydx, 1e-8, x) # why 1e-8 instead of zero?

```


```python
plt.plot(x,y_exact,lw=2, c='r', label='exact+')
plt.plot(x,-y_exact,lw=2, c='r', label='exact-')
plt.plot(x,y_numerical,lw=2, c='k', ls='--',label='numerical+')
plt.plot(x,-y_numerical,lw=2, c='k', ls='--',label='numerical-')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

```

## Example 1: Biochemistry:

The Michaelis-Menten equation relates the rate of an enzyme reaction to its 
substrate concentration, and has the form (if there is a large quantity of 
enzyme):

$${{\rm d}s\over {\rm d}t}=-{V.s\over K+s}$$

where $s=[S]$ is the substrate concentration, and $V$ and $K$ are constants.
By integration, we will derive expressions relating substrate concentration to time when

1. when $s\gg K$
1. when $s\ll K$

and sketch the graphs of $\;s\;$ against $\;t\;$ in each case.

## Case 1

When $s\gg K$, $K$ in the differential equation's denominator can be neglected:

$$\frac{{\rm d}s}{{\rm d}t} \approx -{V.s\over s} = -V$$

The rate of change in concentration is a constant, so this is a **zeroth order process**.

We can integrate directly to obtain
$$ s = \int -V\,{\rm d}t = -Vt + c$$
where $c$ is a constant of integration.

If $s(0)=s_0$, then $c=s_0$.

The graph thus has slope $-V$ and crosses the vertical axis at $s=s_0$.

## Case 2

When $s\ll K$, $\;s\;$ in the differential equation's denominator can be neglected:
$$\frac{{\rm d}s}{{\rm d}t} = -{V.s\over K+s}\quad{\rm~becomes:\quad}\frac{{\rm d}s}{{\rm d}t} = -{V.s\over K}$$

The rate of change in concentration depends on the concentration, so this is a **first order process**.

$$\int {1\over s}\,{\rm d}s = -\int {V\over K}\,{\rm d}t$$

$$\Rightarrow\quad\ln s = -{V\over K}t + D\qquad{\rm D~is~a~constant}$$

$$\Rightarrow\quad s = Be^{-V t\over K}\qquad{\rm B~is~a~constant}$$

On a graph $\;s\;$ crosses the vertical axis at $\;s=B\;$ and decreases with time exponentially.

Lets look at this numerically using the original equation:

$${{\rm d}s\over {\rm d}t}=-{V.s\over K+s}$$


```python
def solve(K):
    V = 1.0
    t0 = 0.0
    t1 = 10.0
    s0 = 1.0
    t = np.linspace(t0,t1,1000)
    def dsdt(s, t):
        return -V * s / (K + s)
    s = scipy.integrate.odeint(dsdt, s0, t)
    plt.plot(t,s)
    plt.xlabel('t')
    plt.ylabel('s')
    
```


```python
interact(solve, K = FloatSlider(value=0.01,min=0.01,max=2.0,step=0.01, continuous_update=False))
```

## Example 2: Bacterial growth:

Suppose the growth rate of a bacterial colony is proportional to its size,
and that we observe that the colony triples in size after 10 hours.
How long will it take to reach 11 times its original size?

Let $\;N\;$ represent the number of bacteria, and suppose there are $\;N_0\;$ at time $\;t=0$.


We are given that when $\;t=10\;$ hours, $\;N=3N_0$, and need to find the time at which $\;N=11N_0\;$.
The equation of growth is:

$$\frac{{\rm d}N}{{\rm d}t} = kN$$

We can solve this to obtain $\;\ln N = kt + c.\;$ Using the $\;t=0\;$ condition, we have that $\;\ln N_0 = c,\;$ so

$$\ln{N\over N_0} = kt$$

We can use the data from $\;t=10\;$ to find $\;k:\;$

$$\ln{3N_0\over N_0} = k 10 \qquad\Rightarrow\qquad k={\ln 3\over 10}=0.11{\rm~hr^{-1}}$$

To find $t$ when $N=11N_0$:

$$\ln{11N_0\over N_0} = 0.11 t \qquad\Rightarrow\qquad t={\ln 11\over 0.11}=21.8{\rm~hours.}$$

## Example 3: Radioactive decay:

The rate at which a sample decays is proportional to the amount left, i.e.

$$\frac{{\rm d}N}{{\rm d}t} = -\lambda N$$

where $\;N\;$ is the mass of radioactive atoms at time $\;t\;$ and $\;\lambda\;$ is called the *decay constant*.

The element radium (atomic mass=226) has a decay constant of $\;13.6 \times 10^{-12}\;$s$^{-1}$.

Solving the differential equation via separation of variables, we obtain the general solution

$$N = Be^{-\lambda t}  ~~~~~~~~(1)$$

Consider an arbitrary initial condition: at $\;t=t_0\;$, $\;N=N_0.\;$ Then:

$$N_0 = Be^{-\lambda t_0} \qquad\Rightarrow\qquad B = N_0 e^{\lambda t_0}$$

Substituting in (1) gives the solution

$$ N = N_0 e^{-\lambda(t-t_0)} ~~~~~~~~~(2)$$

The **half-life**, $\;t_{1\over 2},\;$ is the time taken for $\;N\;$ to reduce by half.

Putting $\;N=N_0/2\;$ and $\;t=t_0+t_{1\over 2}\;$ in (2) we get

$$\ln {N_0\over2N_0} = -\lambda t_{1\over 2}$$

$$t_{1\over 2} = {\ln 2\over\lambda} \approx {0.693\over\lambda}$$

Note that this time is **independent** of the initial value $N_0$,

The half-life for radium is thus $\;t_{1\over 2}={\ln2\over13.6\times10^{-12}}=5.10\times10^{10}\,s,\;$ or about 1600 years.

## Example 4: More biochemistry:

The power to which the concentration of a species is raised in a rate law
is the **order** of the process/reaction with respect to that species.
A reaction with the rate law

$$\frac{{\rm d}X}{{\rm d}t} = k_2~[A]~[B]$$

is first order in $A$ and first order in $B$.

The **overall order** of a process/reaction is the sum of the orders of all the 
components. This rate law is thus second-order overall.

Note that this is **different** from the order of an ODE, which is given by the 
highest derivative.
Both zeroth and first order *processes* are modelled below by *first order 
differential equations*.

## (A) Zeroth order processes:

- rate of change is **independent** of concentration, i.e. the rate of change is proportional to concentration raised to power zero 

$$\frac{{\rm d}A}{{\rm d}t} = k \quad \text{(growth)}$$

$$\frac{{\rm d}A}{{\rm d}t} = -k\quad \text{(decay)}$$

General solutions:
$$A = A_0 + k(t-t_0)$$

$$A = A_0 - k(t-t_0)$$

## (B) First order processes:

- rate of change depends on the concentration of one species, i.e. the rate of change is proportional to concentration raised to first power.

- Half-life is a constant, i.e. is INDEPENDENT of the amount there at the 
  beginning.

$$\frac{{\rm d}A}{{\rm d}t} = kA\quad\text{growth}$$
$$\frac{{\rm d}A}{{\rm d}t} = -kA\quad\text{decay}$$

General solutions:

$$A = A_0 e^{k(t-t_0)}$$
$$A = A_0 e^{-k(t-t_0)}$$


```python
def first_order_process(k):
    t = np.linspace(0,1,1000)
    
    growth = np.exp(k*t)
    decay = growth[-1]*np.exp(-k*t)
    
    plt.plot(t,growth,label='first order growth')
    plt.plot(t,decay,label='first order decay')
    plt.xlabel('x')
    plt.ylabel('A')
    plt.legend()
    
```


```python
interact(first_order_process, k = FloatSlider(value=0.01,min=2,max=10,step=0.01, continuous_update=False))
```

798ccc6047348d7baff3a47c15fd491b6da36830
ipynb
Jupyter Notebook
lectures/lecture-09-differential-eq1.ipynb
SABS-R3/2020-essential-maths
5a53d60f1e8fdc04b7bb097ec15800a89f67a047
Apache-2.0
2021-11-27T12:07:13.000Z
2021-11-27T12:07:13.000Z
lectures/lecture-09-differential-eq1.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
lectures/lecture-09-differential-eq1.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
2020-10-30T17:34:52.000Z
2020-10-30T17:34:52.000Z
Qwen/Qwen-72B
1. YES
2. YES