# Lecture 14

## Systems of Differential Equations I:

### Linear Systems


```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

import cmath

%matplotlib inline

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

global_fig_width = 8
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
##################################################
```

Our aim is to solve systems of equations of the form:

$$\frac{{\rm d}y_i}{{\rm d}t} = f_i(y_1,\ldots,y_n,t),$$

for $\displaystyle i=1,\ldots,n$.


Let us first consider the simplest form: a $2\times2$ **linear** system

  \begin{eqnarray*}
    \frac{{\rm d}x}{{\rm d}t} &=& ax+by,\\
    \frac{{\rm d}y}{{\rm d}t} &=& cx+dy.
  \end{eqnarray*}

## Motivation

- Lotka-Volterra (predator-prey)

- SIR model (epidemics)

- Chemical reaction networks

- Heart cell dynamics

In order to understand these systems, we must first understand coupled linear systems.

## Recap of eigenvalues

$A = \left(\begin{matrix} 1 & 1\\2 & 0  \end{matrix}\right)$


```python
A = sp.Matrix([[1, 1],[2, 0]])
A.eigenvals()

```

## Recap of eigenvectors

$A = \left(\begin{matrix} 1 & 1\\2 & 0  \end{matrix}\right)$


```python
A.eigenvects()

```

## Recap of diagonalisation

Recall that, for a matrix $A$ with **eigenvectors** $\mathbf{v}_1$ and $\mathbf{v}_2$, and **eigenvalues** $\lambda_1$ and $\lambda_2$, we can write a matrix of eigenvectors: $P = \left(\begin{matrix} \mathbf{v}_1 & \mathbf{v}_2  \end{matrix}\right)$. Then:

$$A = P \left(\begin{matrix} \lambda_1 & 0 \\ 0&
      \lambda_2  \end{matrix}\right) P^{-1}$$

(This is also true for general $n \times n$ matrices.)

In our example, $\qquad P = \left(\begin{matrix} -\frac{1}{2} & 1 \\
                                1& 1 \end{matrix}\right),\qquad\lambda_1 = -1,\qquad\lambda_2=2.$

## Coupled ODEs

For coupled system of first order linear differential equations of the form

  \begin{eqnarray*}
    \frac{{\rm d}x}{{\rm d}t} &=& ax+by,\\
    \frac{{\rm d}y}{{\rm d}t} &=& cx+dy.
  \end{eqnarray*}

we have three methods of analysing them mathematically:
  - Turn them into one second order equation (if we can solve second order!)
  - Divide one by other, to get one equation independent of $t$
  - Perform matrix diagonalisation (extends to $n \times n$ problems)

## Example

Solve
  \begin{eqnarray*}
    \frac{{\rm d}x}{{\rm d}t} &=& x+y,\\
    \frac{{\rm d}y}{{\rm d}t} &=& 2x.
  \end{eqnarray*}


Subject to
  \begin{eqnarray*}
    x(0) &=& 0,\\
    y(0) &=& 3.
  \end{eqnarray*}

## Method 1: Second order

We start with:

  \begin{eqnarray*}
    \frac{{\rm d}x}{{\rm d}t} &=& x+y,\\
    \frac{{\rm d}y}{{\rm d}t} &=& 2x.
  \end{eqnarray*}
  
We can convert that into a second order equation:

$$
\frac{{\rm d}^2x}{{\rm d}t^2}
=
\frac{{\rm d}x}{{\rm d}t} + \frac{{\rm d}y}{{\rm d}t}
=
\frac{{\rm d}x}{{\rm d}t} + 2x
\quad \quad\implies \quad\quad
\boxed{\frac{{\rm d}^2x}{{\rm d}t^2}
=
\frac{{\rm d}x}{{\rm d}t} + 2x}
$$

## Method 2: eliminate $t$

We start with:

  \begin{eqnarray*}
    \frac{{\rm d}x}{{\rm d}t} &=& x+y,\\
    \frac{{\rm d}y}{{\rm d}t} &=& 2x.
  \end{eqnarray*}

Then, dividing:

$$\frac{{\rm d}y}{{\rm d}x}
=
\frac{ \quad \frac{{\rm d}y}{{\rm d}t} \quad }{ \frac{{\rm d}x}{{\rm d}t} }
\quad \quad\implies \quad\quad
\boxed{\frac{{\rm d}y}{{\rm d}x}
=
\frac{2x}{x+y} }
$$

## Method 3: diagonalisation

Let $\mathbf{v}=\left(\begin{matrix}x\\y\end{matrix}\right),$


then
$$
    \frac{{\rm d}x}{{\rm d}t} = x+y,\quad
    \frac{{\rm d}y}{{\rm d}t} = 2x, \quad \implies \quad \boxed{\frac{{\rm d}\mathbf{v}}{{\rm d}t} = A\mathbf{v},} 
$$


where  $A = \left(\begin{matrix} 1 & 1\\2 & 0  \end{matrix}\right).$



Substitute
$$
A = P \left(\begin{matrix} \lambda_1 & 0 \\ 0 &
      \lambda_2  \end{matrix}\right) P^{-1}.
$$


then
$$
\frac{{\rm d}\mathbf{v}}{{\rm d}t} = P \left(\begin{matrix} \lambda_1 & 0 \\ 0 &
       \lambda_2  \end{matrix}\right) \left( P^{-1}\mathbf{v} \right)\qquad       
$$


or


$$
\frac{{\rm d}}{{\rm d}t}(P^{-1} \mathbf{v}) = \left(\begin{matrix} \lambda_1 & 0 \\ 0 &
       \lambda_2  \end{matrix}\right) \left( P^{-1}\mathbf{v} \right)
$$


We can now introduce a new variable $\;\mathbf{z} =    \left(\begin{matrix} z_1 \\ z_2\end{matrix}\right)   = P^{-1} \mathbf{v}\;$ so that:


$$
\frac{{\rm d}\mathbf{z}}{{\rm d}t} = \left(\begin{matrix} \lambda_1 & 0 \\ 0 &
       \lambda_2  \end{matrix}\right) \mathbf{z}.
$$


But now, because the matrix is diagonal, the system is no longer coupled.
The first equation **only** involves $z_1$ and the second **only** involves $z_2$, so we can solve each one individually:

$$\frac{{\rm d}z_1}{{\rm d}t} = \lambda_1 z_1 \qquad\implies\qquad z_1(t) = A\,e^{\lambda_1 t}$$


$$\frac{{\rm d}z_2}{{\rm d}t} = \lambda_2 z_2 \qquad\implies\qquad z_2(t) = B\,e^{\lambda_2 t}$$

Finally, we can substitute $z_1$ and $z_2$ back in terms of $x$ and $y$ to find the solution to the original coupled system:


We have

$$\mathbf{z} = \left(\begin{matrix} z_1 \\ z_2\end{matrix}\right) \quad = \quad P^{-1} \mathbf{v} \quad = \quad \frac{1}{3}\left(\begin{matrix} -2 & 2 \\ 2& 1  \end{matrix}\right) \left(\begin{matrix} x \\ y\end{matrix}\right) \quad = \quad \left(\begin{matrix} A\,e^{\lambda_1 t} \\ B\,e^{\lambda_2 t}\end{matrix}\right) $$

Rearranging, we now have two equations relating $x$ and $y$:

$$-2x + 2y = C\,e^{\lambda_1 t}$$

$$2x + y = D\,e^{\lambda_2 t}$$


where $C=3A$ and $D=3B$.  Using our initial conditions, $\,x(0)=0\,$ and $\,y(0)=3\,$ we find $C=6$ and $D=3$.

Finally, solving the simultaneous equations, we have a solution:

$$x(t) = -e^{\lambda_1 t} + e^{\lambda_2 t}$$

$$y(t) = 2e^{\lambda_1 t} + e^{\lambda_2 t}$$

## Summary

- Three methods to analytically solve systems of linear first
  order ODEs
- Best method depends on the system and what you need to ask about it


eda22f6708393ebc59220d1882e037605fb8902e
ipynb
Jupyter Notebook
lectures/lecture-14-systems1.ipynb
SABS-R3/2020-essential-maths
5a53d60f1e8fdc04b7bb097ec15800a89f67a047
Apache-2.0
2021-11-27T12:07:13.000Z
2021-11-27T12:07:13.000Z
lectures/lecture-14-systems1.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
lectures/lecture-14-systems1.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