# 💃 Classifying differential equations 💃



Questions:
- What is a <mark>differential equation</mark>?
- What is the difference between an ordinary (<mark>ODE</mark>) and partial (<mark>PDE</mark>) differential equation? 
- How do I classify the different types of differential equations?

Objectives:
- Identify the dependent and independent variables in a differential equation
- Distinguish between and ODE and PDE
- Identify the <mark>order</mark> of a differential equation
- Distinguish between <mark>linear</mark> and <mark>non-linear</mark> equations
- Distinguish between <mark>heterogeneous</mark> and <mark>homogeneous</mark> equations
- Identify a <mark>separable</mark> equation

### A differential equation is an equation that relates one or more functions and their derivatives

- The functions usually represent physical quantities (e.g. $\mathbf{F}$))
- The derivative represents a rate of change (e.g. acceleration)
- The differential equation represents the relationship between the two.
- For example, Newton's second law for $n$ particles of mass $m$:

\begin{equation}
\mathbf{F}(t,\mathbf{x},\mathbf{v}) = m\frac{d\mathbf{v}}{dt}
\end{equation}

### An independent variable is... a quantity that varies independently...

- An <mark>independent variable</mark> does not depend on other variables
- A <mark> dependent variable</mark> depends on the independent variable

\begin{equation}
\mathbf{F}(t,\mathbf{x},\mathbf{v}) = m\frac{d\mathbf{v}}{dt}
\end{equation}

- $t$ is an <mark>independent</mark> variable
- $x$ and $v$ are <mark>dependent</mark> variables
- Writing $x = x(t)$ makes this relationship clear.
 

### Differential equations can be classified in a variety of ways

There are several ways to describe and classify differential equations. There are standard solution methods for each type, so it is useful to understand the classifications. 



<small>Once you can cook a single piece of spaghetti, you can cook all pieces of spaghetti!</small>

### An ODE contains differentials with respect to only one variable

For example, the following equations are ODEs:

\begin{equation}
\frac{d x}{d t} = at
\end{equation}
\begin{equation}
\frac{d^3 x}{d t^3} + \frac{x}{t} = b
\end{equation}

As in each case the differentials are with respect to the single variable $t$.



### Partial differential equations (PDE) contain differentials with respect to several independent variables. 

An example of a PDE is:

\begin{equation}
\frac{\partial x}{\partial t} = \frac{\partial x}{\partial y} 
\end{equation}

As there is one differential with respect to $t$ and one differential with respect to $y$. 

Note also the difference in notation - <mark>for ODEs we use $d$ whilst for PDEs we  use $\partial$</mark>.



### The order of a differential equation is the highest order of any differential contained in it. 

For example:

$\frac{d x}{d t} = at$ is <mark>first order</mark>. 

$\frac{d^3 x}{d t^3} + \frac{x}{t} = b$ is <mark>third order</mark>.

**Note:** $\frac{d^3 x}{d t^3}$ does not equal $\left(\frac{d x}{d t}\right)^3$!

###  Linear equations do not contain higher powers of either the dependent variable or its differentials

For example:

$\frac{d^3 x}{d t^3} = at$ and $\frac{\partial x}{\partial t} = \frac{\partial x}{\partial y} $ are <mark>linear</mark>.

$(\frac{d x}{d t})^3 = at$ and $\frac{d^3 x}{d t^3} = x^2$ are <mark>non-linear</mark>.

Non-linear equations can be particularly nasty to solve analytically, and so are often tackled numerically.




### Homogeneous equations do not contain any non-differential terms

For example:

$\frac{\partial x}{\partial t} = \frac{\partial x}{\partial y}$ is a <mark>homogeneous equation</mark>.

$\frac{\partial x}{\partial t} - \frac{\partial x}{\partial y}=a$ is a <mark>heterogeneous equation</mark> (unless $a=0$!).



### Separable equations can be written as a product of two functions of different variables

A separable differential equation takes the form 

\begin{equation}
f(x)\frac{d x}{d t} = g(t)
\end{equation}

Separable equations are some of the easiest to solve as we can split the equation into two independent parts with fewer variables, and solve each in turn - we will see an example of this in the next tutorial.


Keypoints:

- An independent variable is a quantity that varies independently
- Differential equations can be classified in a variety of ways
- An ODE contains differentials with respect to only one variable
- The order is the highest order of any differential contained in it
- Linear equations do not contain higher powers of either the dependent variable or its differentials
- Homogeneous equations do not contain any non-differential terms

# Radioactive Decay


----
Questions:
- How can I describe radioactive decay using a first-order ODE? 
- What are <mark>initial conditions</mark> and why are  they important?
----

----
Objectives:
- Map between physical notation for a particular problem and the more general notation for all differential equations
- Solve a linear, first order, separable ODE using integration
- Understand the physical importance of initial conditions 
----



### Radioactive decay can be modelled a linear, first-order ODE

As our first example of an ODE we will model radioactive decay using a differential equation.

We know that the decay rate is proportional to the number of atoms present.  Mathematically, this relationship can be expressed as:

\begin{equation}
\frac{d N}{d t} = -\lambda N
\end{equation}

Note that we could choose to use different variables, for example:

\begin{equation}
\frac{d y}{d x} = cy
\end{equation}

However we try to use variables connected to the context of the problem. For example $N$ for the Number of atoms.


If we know that 10% of atoms will decay per second we could write:

\begin{equation}
\frac{d N}{d t} = -0.1 N
\end{equation}

where $N$ is the number of atoms and $t$ is time measured in seconds.

This equation is linear and first-order.

| physical notation | generic notation $\left(\frac{dy}{dx} = f(y)\right)$ |
|-----|-----|
|number of  atoms $N$ | dependent  variable $y$|
| time $t$ | independent variable  $x$|
| decay rate $\frac{dN}{dt}$ | differential $\frac{dy}{dx}$|
| constant of proportionality $\lambda=0.1 $ | parameter(s) |

### The equation for radioactive decay is separable and has an analytic solution

The radioactive decay equation is <mark>separable</mark>:

\begin{equation}
\frac{1}{N}\frac{d N}{d t} = -\lambda 
\end{equation}

\begin{equation}
\frac{dN}{N} = -\lambda dt.
\end{equation}

We can then integrate each side:

\begin{equation}
\ln N = -\lambda t + const.
\end{equation}

and solve for N:

\begin{equation}
N = e^{-\lambda t}e^{\textrm{const.}} 
\end{equation}


### To model a physical system an initial value has to be provided

\begin{equation}
N(t) = e^{-\lambda t}e^{\textrm{const.}} 
\end{equation}

At the beginning (when $t=0$):

\begin{equation}
N (t_0) = e^{0}e^{\textrm{const.}} = e^{\textrm{const.}}
\end{equation}

So we can identify $e^{\textrm{const.}}$ as the amount of radioactive material that was present in the beginning. We denote this starting amount as $N_0$.

Substituting this back into Equation 4, the final solution can be more meaningfully written as:

\begin{equation}
N = N_0 e^{-\lambda t}
\end{equation}

We now have  not  just one solution, but a whole class  of solutions that are dependent on the initial amount of radioactive material $N_0$. 

Remember that not all mathematical solutions make physical sense. <mark>To model a physical system, this initial value (also known as initial condition) has to be provided alongside the constant of proportionality $\lambda$.</mark>



### ODEs can have initial values or boundary values

ODEs have either initial values or boundary values. Take Newton's second law as an example:

\begin{equation}
\frac{\mathrm{d}^2x}{\mathrm{d}t^2} = -g
\end{equation}

An initial value problem would be where we know the starting position and velocity. A boundary value problem would be where we specify the position of the ball at times $t=t_0$ and $t=t_1$. 

For second-order ODEs (such as acceleration under gravity) we need to provide two initial/boundary conditions, for third-order ODEs we would need to provide three, and so on. Our radioactive decay example was a first-order ODE and so we only had to provide a single initial condition.


-----
Keypoints:
- Radioactive decay can be modelled a linear, first-order ODE
- The equation for radioactive decay is separable and has an analytic solution
- To model a physical system an initial value has to be provided
- The number of  initial conditions depends on the order of  the differential equation
-----

# Euler's Method

-------
Questions:
- How do I use Euler's method to solve a first-order ODE?



--------



---------
Objectives:
- Use Euler's method, implemented in Python, to solve a first-order ODE
- Understand that this method is approximate and the significance of step size $h$
- Compare results at different levels of approximation using the `matplotlib` library.
------

### There are a variety of ways to solve an ODE

In the previous lesson we considered nuclear decay:

\begin{equation}
\frac{\mathrm{d} N}{\mathrm{d} t} = -\lambda N
\end{equation}

This is one of the simplest examples of am ODE - a first-order, linear, separable differential equation with one dependent variable. We saw that we could model the number of atoms $N$ by finding an analytic solution through integration:

\begin{equation}
N = N_0 e^{-\lambda t}
\end{equation}

However there is more than one way to crack an egg (or solve a differential equation). We could have, instead, used an approximate, numerical method. One such method - Euler's method - is this subject of this lesson.

### A function can be approximated using a Taylor expansion

The Taylor series is a polynomial expansion of a function about a point.
For example, the image below shows $\mathrm{sin}(x)$ and its Taylor approximation by polynomials of degree 1, 3, 5, 7, 9, 11, and 13 at $x = 0$. 
<figure>

</figure>
Credit: Image By IkamusumeFan - CC BY-SA 3.0, Commons Wikimedia


The Taylor series of $f(x)$ evaluated at point $a$ can be expressed as:

\begin{equation}
f(x) = f(a) + \frac{\mathrm{d} f}{\mathrm{d} x}(x-a) + \frac{1}{2!} \frac{\mathrm{d} ^2f}{\mathrm{d} x^2}(x-a)^2 + \frac{1}{3!} \frac{\mathrm{d} ^3f}{\mathrm{d} x^3}(x-a)^3
\end{equation}

Returning to our example of nuclear decay, we can use a <mark>Taylor expansion</mark> to write the value of $N$ a short interval $h$ later:

\begin{equation}
N(t+h) = N(t) + h\frac{\mathrm{d}N}{\mathrm{d}t} + \frac{1}{2}h^2\frac{\mathrm{d}^2N}{\mathrm{d}t^2} + \ldots
\end{equation}

\begin{equation}
N(t+h) = N(t) + hf(N,t) + \mathcal{O}(h^2)
\end{equation}

If $h$ is small and $h^2$ is very small we can neglect the terms in $h^2$ and higher and we get:

\begin{equation}
N(t+h) = N(t) + hf(N,t).
\end{equation}


```python
# https://gist.github.com/dm-wyncode/55823165c104717ca49863fc526d1354
"""Embed a YouTube video via its embed url into a notebook."""
from functools import partial

from IPython.display import display, IFrame

width, height = (560, 315, )



def _iframe_attrs(embed_url):
    """Get IFrame args."""
    return (
        ('src', 'width', 'height'), 
        (embed_url, width, height, ),
    )



def _get_args(embed_url):
    """Get args for type to create a class."""
    iframe = dict(zip(*_iframe_attrs(embed_url)))
    attrs = {
        'display': partial(display, IFrame(**iframe)),
    }
    return ('YouTubeVideo', (object, ), attrs, )

def youtube_video(embed_url):
    """Embed YouTube video into a notebook.

    Place this module into the same directory as the notebook.

    >>> from embed import youtube_video
    >>> youtube_video(url).display()
    """
    YouTubeVideo = type(*_get_args(embed_url)) # make a class
    return YouTubeVideo() # return an object

```


```python

youtube_video("https://www.youtube.com/embed/3d6DsjIBzJ4").display()
```







### Euler's method can be used to approximate the solution of differential equations

We can keep applying the equation above so that we calculate $N(t)$ at a succession of equally spaced points for as  long as we want. If $h$ is small enough we can get a good approximation to the solution of the equation. This method for solving differential equations is called Euler's method, after Leonhard Euler, its inventor.



**Note**: Although we are neglecting terms $h^2$ and higher, Euler's method typically has an error linear in $h$ as the error accumulates over repeated steps. This means that if we want to double the accuracy of our calculation we need to double the number of steps, and double the calcuation time.


### Euler's method can be applied using the Python skills we have developed

Let's use Euler's method to solve the differential equation for nuclear decay. We will model the decay process over a period of 10 seconds, with the decay constant $\lambda=0.1$ and the initial condition $N_0 = 1000$.

\begin{equation}
\frac{\mathrm{d}N}{\mathrm{d} t} = -0.1 N
\end{equation}


----

Keypoints:

- There are a variety of ways to solve an ODE
- A function can be approximated using a Taylor expansion
- If the step size $h$ is small then higher  order terms can be neglected
- Euler's method can be used to approximate the solution of differential equations
- Euler's method can be applied using the Python skills we have developed
- We can easily visualise our results, and compare against the analytical solution, using the matplotlib plotting library
-----

# The Strange Attractor 

Questions:
- How do I solve simultaneous ODEs?
- How do I solve second-order ODEs (and higher)?


Objectives:
- Use Euler's method, implemented in Python, to solve a set of simultaneous ODEs
- Use Euler's method, implemented in Python, to solve a second-order ODE
- Understand how the same method could be applied to higher order ODEs

### Computers don't care so much about the type of differential equation

- In the previous lesson we used Euler's method to model radioactive decay. But we have seen that this equation can also be solved analytically

- However there are a large number of physical equations that cannot be solved analytically, and that rely on numerical methods for their modelling. 

- For example, the Lotka-Volterra equations for studying predator-prey interactions have multiple dependent variables and the Cahn-Hilliard equation for modelling phase separation in fluids is non-linear. These equations can be solved analytically for particular, special cases only.

- Computers don't care whether an equation is linear or non-linear, multi-variable or single-variable, the numerical method for studying it is much the same 

- The caveat is that numerical methods are approximate and so we need to think about the accuracy of the methods used.



### Simultaneous ODEs can be solved using numerical methods

One class of problems that are difficult to solve analytically are <bold>simultaenous ODEs</bold>. These are equations where the derivative of each dependent variable can depend on any of the variables (dependent or independent). For example,

\begin{eqnarray}
\frac{\mathrm{d}x}{\mathrm{d}t} &=& xt + y \\
\frac{\mathrm{d}y}{\mathrm{d}t} &=& \mathrm{sin}^2\omega t - xy
\end{eqnarray}


### Euler's method is easily extended to simultaneous ODEs

In the previous lesson we were introduced to Euler's method for the single variable case:

\begin{equation}
N(t+h) = N(t) + hf(N,t).
\end{equation}

This can be easily extended to the multi-variable case using vector notation:

\begin{equation}
\mathbf{r}(t+h) = \mathbf{r}(t) + h\mathbf{f}(\mathbf{r},t).
\end{equation}

We have seen that arrays can be easily represented in Python using the NumPy library. This allows us to do arithmetic with vectors directly (rather than using verbose workarounds such as `for` loops), so the code is not much more complicated than the one-variable case.

### The Lorenz equations are a famous set of simultaneous ODEs

The Lorenz equations demonstrate that systems can be deterministic and yet be inherently unpredictable (even in the absence of quantum effects).

\begin{eqnarray}
\frac{\mathrm{d}x}{\mathrm{d}t} &=& \sigma(y-x) \\
\frac{\mathrm{d}y}{\mathrm{d}t} &=& rx-y-xz \\
\frac{\mathrm{d}z}{\mathrm{d}t} &=& xy-bz
\end{eqnarray}

There three dependent variables - $x$, $y$ and $z$, and one independent variable $t$. There are also three constants -  $\sigma$, $r$ and $b$.

For particular values of $\sigma$, $r$ and $b$ the Lorenz systems has <bold>chaotic behaviour</bold> (a strong dependence on the initial conditions) and for a subset of these, there is also fractal structure called the <bold>strange attractor</bold>.

### Higher order ODEs can be re-cast as simultaneous ODEs and solved in the same way

Many physical equations are second-order or higher. The general form for a second-order differential equation with one dependent variable is:

\begin{equation}
\frac{\mathrm{d}^2x}{\mathrm{d}t^2} = f\left(x,\frac{\mathrm{d}x}{\mathrm{d}t},t\right)
\end{equation}

Luckily, we can re-cast a higher order equation as a set of simultaneous equations, and then solve in the same way as above.



Let's use the non-linear pendulum as an example. For a pendulum with an art of length $l$ and a bob of mass $m$, Newton's second law ($F=ma$) provides the following equation of motion:

\begin{equation}
ml\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = -mg\sin(\theta),
\end{equation}

where $\theta$ is the angle of displacement of the arm from the vertical and $l\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2}$ is the acceleration of the mass in the tangential direction. The exact solution to this equation is unknown, but we now have the knowledge needed to find a numerical approximation.



The equation can be re-written as:

\begin{equation}
\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = -\frac{g}{l}\sin(\theta)
\end{equation}

We define a new variable $\omega$:

\begin{equation}
\frac{\mathrm{d}\theta}{\mathrm{d}t} = \omega
\end{equation}

and substitute this into the equation of motion:

\begin{equation}
\frac{\mathrm{d}\omega}{\mathrm{d}t} = -\frac{g}{l}\sin(\theta).
\end{equation}

These two expressions form a set of simultaneous equations that can be solved numerically using the method outlined above.


-----
Keypoints:
- Computers don't care so much about the type of differential equation
- Simulatenous ODEs can also be solved using numerical methods
- Euler's method is easily extended to the multi-variable case
- Higher order ODEs can be re-cast as simultaneous ODEs and solved the same way
-----

# The Runge-Kutta method




Questions:
- How do I use the Runge-Kutta method for more accurate solutions?




Objectives:
- Use the Runge-Kutta method, implemented in Python, to solve a first-order ODE
- Compare results at different levels of approximation using the `matplotlib` library.

### The Runge-Kutta method is more accurate than Euler's method and runs just as fast

So far we have used Euler's method for solving ODEs. We have learnt that, using this method, the final expression for the total error is linear in $h$. 

However for roughly the same compute time we can reduce the total error so it is of order $h^2$ by implementing another method - the <bold>Runge-Kutta method</bold>.



> Note: It is common to use the Runge-Kutta method for solving ODEs given the improved accuracy over Euler's method. However Euler's method is still commonly used for PDEs (where there are other, larger, sources of error).


> Note: The Runge-Kutta method is actually a family of methods. In fact, Euler's method is the first-order Runge-Kutta method. There is then the second-order Runge-Kutta method, third-order Runge-Kutta method, and so on..

Euler's method does not take into account the curvature of the solution, whilst Runge-Kutta methods do, by calculating the gradient at intermediate points in the (time-)step. 



### Runge-Kutta methods are derived from Taylor expansion(s) around intermediate point(s)

To derive the second-order Runge-Kutta method we:

1) estimate $x(t+h)$ using a Taylor expansion around $t+\frac{h}{2}$:

\begin{equation}
x(t+h) = x(t+\frac{h}{2}) + \frac{h}{2}\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right)_{t+\frac{h}{2}} + \frac{h^2}{8}\left(\frac{\mathrm{d}^2x}{\mathrm{d}t^2}\right)_{t+\frac{h}{2}}+\mathcal{O}(h^3)
\end{equation}

2) estimate $x(t)$ using a Taylor expansion around $t-\frac{h}{2}$:

\begin{equation}
x(t) = x(t+\frac{h}{2}) - \frac{h}{2}\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right)_{t+\frac{h}{2}} + \frac{h^2}{8}\left(\frac{\mathrm{d}^2x}{\mathrm{d}t^2}\right)_{t+\frac{h}{2}}+\mathcal{O}(h^3)
\end{equation}

3) Subtract Equation 2 from Equation 1 and re-arrange:

\begin{eqnarray}
x(t+h) &=& x(t) + h\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right))_{t+\frac{h}{2}}+\mathcal{O}(h^3) \\
x(t+h) &=& x(t) + hf(x(t+\frac{h}{2}),t+\frac{h}{2})+\mathcal{O}(h^3)
\end{eqnarray}



> Note: the $h^2$ term has completely disappeared, and the error term is now order $h^3$. We can say that this approximation is now accurate to order $h^2$.



The problem is that this requires knowledge of $x(t+\frac{h}{2})$ which we don't currently have. We can however estimate this using the Euler method!

\begin{equation}
x(t+\frac{h}{2}) = x(t) + \frac{h}{2}f(x,t).
\end{equation}

Substituting this into Equation 3 above, we can write the method for a single step as follows:

\begin{eqnarray}
k_1 &=& hf(x,t) \\
k_2 &=& hf(x+\frac{k_1}{2},t+\frac{h}{2})\\
x(t+h) &=& x(t) + k_2
\end{eqnarray}

See how $k_1$ is used to give an estimate for $x(t+\frac{h}{2})$ ($k_2$), which is then substituted into the third equation to give an estimate for  $x(t+h)$.

> Note: Higher orde Runge-Kutta methods can be derived in a similar way - by calculating the Taylor series around various points and then taking a linear combination of the resulting expansions. As we increase the number of intermediate points, we increase the accuracy of the method. The downside is that the equations get increasingly complicated. 


Keypoints:
- The Runge-Kutta method is more accurate than Euler's method and runs just as fast
- Runge-Kutta methods are derived from Taylor expansion(s) around intermediate point(s)
- Runge-Kutta methods can be applied using the Python skills we have developed
- We can easily compare our various models using the `matplotlib` plotting library

53713b68fad0443e7a00b74d5e8b84b87bc7ee6b
ipynb
Jupyter Notebook
slides/ODE-slides.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
slides/ODE-slides.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-10-06T08:11:43.000Z
2021-11-24T13:46:08.000Z
slides/ODE-slides.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-11-03T10:11:28.000Z
2021-11-03T10:11:28.000Z
Qwen/Qwen-72B
1. YES
2. YES