![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "fLE4p9FNAuja"
},
"source": [
"# Differential Equations"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "65f7lpUoAujd"
},
"source": [
"Differential equations is without doubt one of the more useful branch of mathematics in science. They are used to model problems involving the change of some variable with respect to another. Differential equations cover a wide range of different applications, ranging from ordinary differential equations (ODE) until boundary-value problems involving many variables. For the sake of simplicity, throughout this section we shall cover only ODE systems as they are more elemental and equally useful. First, we shall cover first order methods, then second order methods and finally, system of differential equations.\n",
"\n",
"See: \n",
"* http://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture21.pdf\n",
"* http://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture22.pdf\n",
"* Infectious Disease Modelling:\n",
" * [Understanding the models that are used to model Coronavirus](https://towardsdatascience.com/infectious-disease-modelling-part-i-understanding-sir-28d60e29fdfc): Explaining the background and deriving the formulas of the SIR model from scratch. Coding and visualizing the model in Python.\n",
" * [Beyond the Basic SIR Model](https://towardsdatascience.com/infectious-disease-modelling-beyond-the-basic-sir-model-216369c584c4)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "_UWnGRIfAujf"
},
"source": [
"- - -\n",
"- [First Order Methods](#First-Order-Methods) \n",
" - [Euler's method](#Euler's-method)\n",
" - [Example 1](#Example-1)\n",
"- [High Order Methods](#High-Order-Methods)\n",
" - [Second-order Runge-Kutta methods](#Second-order-Runge-Kutta-methods)\n",
" - [Example 2](#Example-2)\n",
" - [Fourth-order Runge-Kutta method](#Fourth-order-Runge-Kutta-method)\n",
"- [Two-Point Boundary Value Problems](#Two-Point-Boundary-Value-Problems)\n",
" - [Example 3](#Example-3)\n",
"- [SIR Model](./material/)\n",
"- - -"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "gECD72StAujg",
"outputId": "2936d89e-4cc2-4f07-f63d-98fdb39d4e57"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "Gs3xyp5_Aujp"
},
"outputs": [],
"source": [
"import numpy as np\n",
"## JSAnimation import available at https://github.com/jakevdp/JSAnimation\n",
"#from JSAnimation import IPython_display\n",
"from matplotlib import animation"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "rXIj_PTAAujt"
},
"source": [
"- - - "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ng_1qE8bAujv",
"tags": []
},
"source": [
"## Physical motivation\n",
"\n",
"### Momentum\n",
"In absence a forces a body moves freely to a constant velocity (See figure)\n",
"\n",
"The quantity which can be associated with the change of speed of a body is the instantaneous _momentum_ of a particle of mass $m$ moving with instantaneous velocity $\\boldsymbol{v}$, given by\n",
"\\begin{align}\n",
" \\boldsymbol{p}=&\\gamma m \\mathbf{v}\\,, &\\text{where: } \\gamma=&\\frac{1}{\\sqrt{1-{|\\boldsymbol{v}|^2}/{c^2}}}\\,,\n",
"\\end{align}\n",
"and $c\\approx 3\\times 10^8\\ $m/s is the speed of light.\n",
"\n",
"If $\\gamma\\approx1$, or equivalent, sif $|\\boldsymbol{v}|\\ll c$:\n",
"\\begin{align}\n",
" \\boldsymbol{p}\\approx &m \\boldsymbol{v}\n",
"\\end{align}\n",
"### Equation of motion\n",
"Any change of the momentum can be atribuited to some force, $\\boldsymbol{F}$. In fact, the Newton's second law can be defined in a general way as\n",
"> The change of the momentum of particle is equal to the net force acting upon the system times the duration of the interaction\n",
"\n",
"For one instaneous duration $\\Delta t$, this law can be written as\n",
"$$\n",
"\\Delta \\boldsymbol{p} = \\boldsymbol{F}\\Delta t,\\qquad \\Delta t \\to 0 \n",
"$$\n",
"or as the _equation of motion_, which is the _differential equation_\n",
"\\begin{align}\n",
"\\boldsymbol{F}=&\\lim_{t\\to 0} \\frac{\\Delta \\boldsymbol{p}}{\\Delta t} \\\\\n",
"\\boldsymbol{F}=&\\frac{\\operatorname{d}\\boldsymbol{p}}{\\operatorname{d} t} \\,.\n",
"\\end{align}\n",
"This equation of motion is of general validity and can be applied numerically to solve any problem.\n",
"\n",
"### Constant mass\n",
"In the special case of constant mass\n",
"\\begin{align}\n",
"\\boldsymbol{F}=&m \\boldsymbol{a}\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Example: Drag force\n",
"For a body falling in the air, in addition to the gravitiy force $\\boldsymbol{W}$, there is a dragging force in the opposite direction given by [[see here for details]](https://en.wikipedia.org/wiki/Drag_(physics))\n",
"$$\n",
"F_D=\\frac{1}{2}{\\rho A C_d} v^2(t)=\\frac{1}{2m^2}{\\rho A C_d} p^2(t)\\,,\n",
"$$\n",
"where \n",
"* $A$: is the frontal area on a plane perpendicular to the direction of motion—e.g. for objects with a simple shape, such as a sphere\n",
"* $\\rho$ is the density of the air\n",
"* $C_d$: _Drag coefficient_. $C_d=0.346$ for a [baseball](http://baseball.physics.illinois.edu/DragTPTMay2014.pdf)\n",
"* $v(t)$: speed of the baseball\n",
"* $p(t)$: momentum of baseball\n",
"\n",
"The equation of motion is\n",
"\\begin{align}\n",
"W-F_D=& \\frac{d p}{dt}\\\\\n",
"\\frac{d p}{dt}=&m g - \\frac{\\rho A C_d}{2m^2} p^2(t). \n",
"\\end{align}\n",
"\n",
"We see then that, in general\n",
"$$\n",
"\\frac{d p}{dt}=f(t,p)\\,.\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "NasuV5ArAujx"
},
"source": [
"\n",
"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Z6DwVvzgAujy"
},
"source": [
"## Mathematical background\n",
"Ordinary differential equations normally implies the solution of an initial-value problem, i.e., the solution has to satisfy the differential equation together with some initial condition. Real-life problems usually imply very complicated problems and even non-soluble ones, making any analytical approximation unfeasible. Fortunately, there are two ways to handle this. First, for almost every situation, it is generally posible to simplify the original problem and obtain a simpler one that can be easily solved. Then, using perturbation theory, we can perturbate this solution in order to approximate the real one. This approach is useful, however, it depends very much on the specific problem and a systematic study is rather complicated.\n",
"\n",
"The second approximation, and the one used here, consists of a complete numerical reduction of the problem, solving it exactly within the accuracy allowed by implicit errors of the methods. For this part, we are going to assume well-defined problems, where solutions are expected to be well-behaved."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "iF1Wi1W6Aujz",
"tags": []
},
"source": [
"## Euler's method"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"### First order differential equations"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "nUie936CAuj0"
},
"source": [
"This method is the most basic of all methods, however, it is useful to understand basic concepts and definitions of differential equations problems. \n",
"\n",
"Suppose we have a well-posed initial-value problem given by:\n",
"\n",
"$$ \\frac{dy}{dt}=f(t,y),\\ \\ \\ a\\leq t\\leq b, \\ \\ \\ \\ y(a) = \\alpha $$\n",
"\n",
"Now, let's define a **mesh points** as a set of values $\\{t_i\\}$ where we are going to approximate the solution $y(t)$. These points can be equal-spaced such that\n",
"\n",
"$$ t_i = a+ i \\Delta t,\\ \\ \\ \\ \\mbox{with}\\ \\ i=1,2,3,\\cdots,N \\ \\ \\mbox{and}\\ \\Delta t = \\frac{b-a}{N}. $$\n",
"\n",
"Here, $h$ is called the **step size** of the **mesh points**.\n",
"\n",
"Now, using the first-order approximation of the derivative that we saw in [Numerical Differentiation](http://nbviewer.ipython.org/github/sbustamante/ComputationalMethods/blob/master/material/numerical-calculus.ipynb#Numerical-Differentiation), we obtain\n",
"\n",
"$$ \\frac{dy}{dt}\\approx \\frac{y(t+\\Delta t)-y(t)}{\\Delta t} $$\n",
"\n",
"or\n",
"\n",
"$$ \\left.\\frac{dy}{dt}\\right|_{t=t_i}\\approx \\frac{y(t_{i+1})-y(t_i)}{\\Delta t} $$\n",
"\n",
"The original problem can be rewritten as\n",
"\n",
"$$ y_{i+1} = y_i + f(t_i, y_i)\\Delta t \\pm \\frac{\\Delta t^2}{2}y^{''}(\\xi_i) $$\n",
"\n",
"where the notation $y(t_i)\\equiv y_i$ has been introduced and the last term (error term) can be obtained taking a second-order approximation of the derivative. $\\xi_i$ is the value that give the maximum in absolute value for $y''$\n",
"\n",
"This equation can be solved recursively as we know the initial value $y_0 = y(a) = \\alpha$."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "9FEdijqEAuj1"
},
"source": [
"#### Second order differential equations"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ccehvb8JAuj3"
},
"source": [
"The formalism of the Euler's method can be applied for any system of the form:\n",
"\n",
"$$ \\frac{dy}{dt}=f(t,y),\\ \\ \\ a\\leq t\\leq b, \\ \\ \\ \\ y(a) = \\alpha $$\n",
"\n",
"However, it is possible to extend this to second-order systems, i.e., systems involving second derivatives. Let's suppose a general system of the form:\n",
"\n",
"$$ \\frac{d^2y}{dt^2}+ g(t,y)\\frac{dy}{dt}=f(t,y),\\ \\ \\ a\\leq t\\leq b, \\ \\ \\ \\ y(a) = \\alpha\\ \\ \\mbox{and}\\ y'(a) = \\beta$$\n",
"\n",
"For this system we have to provide both, the initial value $y(a)$ and the initial derivative $y'(a)$.\n",
"\n",
"Now, let's define a new variable $w(t) = y'(t)$, the previous problem can be then written as \n",
"\n",
"\\begin{align}\n",
"\\frac{dy}{dt}=&w(t) \\\\ \n",
"\\frac{dw}{dt}=&f(t,y)-g(t,y)w\\, ,\\ \\ \\ a\\leq t\\leq b\\,, \\ \\ \\ \\ y(a) = \\alpha\\ \\ \\mbox{and}\\ w(a) = \\beta\n",
"\\end{align}\n",
"\n",
"Each of this problem has now the form required by the Euler's algorithm, and the solution can be computed as:\n",
"\n",
"\\begin{align}\n",
"w_{i+1}=& w_{i} + [f(t_i,y_i)-g(t_i,y_i)w_i]\\Delta t \\\\ \n",
"y_{i+1}=&y_i + w_i \\Delta t\\, ,\\ \\ \\ a\\leq t\\leq b\\,, \\ \\ \\ \\ y(a) = \\alpha\\ \\ \\mbox{and}\\ w(a) = \\beta\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Euler method for second order differential equations\n",
"We can define the column vectors\n",
"\\begin{align}\n",
"\\boldsymbol{U}=\\begin{bmatrix}\n",
"U_1\\\\\n",
"U_2\\\\\n",
"\\end{bmatrix}=&\\begin{bmatrix}\n",
"y(t)\\\\\n",
"w(t)\\\\\n",
"\\end{bmatrix}\\,,& \n",
"\\boldsymbol{V}(\\boldsymbol{U},t)=\\begin{bmatrix}\n",
"w(t)\\\\\n",
"f(t,y)-g(t,y)w(t)\\\\\n",
"\\end{bmatrix}=&\\begin{bmatrix}\n",
"U_2\\\\\n",
"f(t,U_1)-g(t,U_1)U_2\\\\\n",
"\\end{bmatrix}\n",
"\\end{align}\n",
"such that\n",
"\n",
"$$\n",
"\\frac{dU}{dt}=V(t,U)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"### Newton's second law of motion"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "2Gru5svcAuj4",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"In fact, the equation of motion can be seen as the system of equations\n",
"\\begin{align}\n",
"\\frac{\\boldsymbol{p}}{m}=&\\frac{\\operatorname{d}\\boldsymbol{x}}{\\operatorname{d}t}\\\\\n",
"\\boldsymbol{p}=& \\gamma m \\boldsymbol{v}\\\\\n",
"\\boldsymbol{F}=&\\frac{\\operatorname{d}\\boldsymbol{p}}{\\operatorname{d} t} \n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "vJRY_2HyAuj5"
},
"source": [
"### Example:\n",
"\n",
"\n",
"An object of 0.5 Kg is launched from the top of a building of 50 m with an horizontal speed of 5 m/s. Study the evolution of the movement\n",
"* Neglecting the air friction"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0. , 0.01])"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"Δt=0.01 #s\n",
"np.arange( 0,3+Δt,Δt)[:2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"\\frac{\\boldsymbol{p}}{m}=&\\frac{\\operatorname{d}\\boldsymbol{x}}{\\operatorname{d}t}\\\\\n",
"\\boldsymbol{p}=& \\gamma m \\boldsymbol{v}\\\\\n",
"\\boldsymbol{F}=&\\frac{\\operatorname{d}\\boldsymbol{p}}{\\operatorname{d} t} \n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "AM__qjSaAuj7"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"df=pd.DataFrame()\n",
"m=.5 ## Kg\n",
"g=9.8 ## m/s^2\n",
"## intial conditions\n",
"x=np.array([0,50,0]) #m\n",
"v=np.array([5,0,0]) ## m/s\n",
"p=m*v\n",
"Δt=0.01 #s\n",
"## Analysis for the first 3 s\n",
"df=df.append({'x':x,'p':p},ignore_index=True)\n",
"Fg=np.array([0,-m*g,0]) #N (Weight)\n",
"\n",
"for t in np.arange( 0,3+Δt,Δt): #s\n",
" p=p+Fg*Δt\n",
" x=x+(p/m)*Δt\n",
" df=df.append({'x':x,'p':p},ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 1882
},
"colab_type": "code",
"id": "4uezYasEAuj_",
"outputId": "165fb893-f4a7-40b1-9301-145c4cb8e496"
},
"outputs": [
{
"data": {
"text/html": [
"| \n", " | x | \n", "p | \n", "
|---|---|---|
| 0 | \n", "[0, 50, 0] | \n", "[2.5, 0.0, 0.0] | \n", "
| 1 | \n", "[0.05, 49.99902, 0.0] | \n", "[2.5, -0.049, 0.0] | \n", "
| 2 | \n", "[0.1, 49.997060000000005, 0.0] | \n", "[2.5, -0.098, 0.0] | \n", "
| \n", " | x | \n", "p | \n", "t | \n", "c | \n", "
|---|---|---|---|---|
| 0 | \n", "[0, 50, 0] | \n", "[2.5, 0.0, 0.0] | \n", "0.00 | \n", "0.0 | \n", "
| 1 | \n", "[0.05, 49.99902, 0.0] | \n", "[2.5, -0.049, 0.0] | \n", "0.01 | \n", "0.0 | \n", "
| 2 | \n", "[0.1, 49.997060000000005, 0.0] | \n", "[2.5, -0.098, 0.0] | \n", "0.02 | \n", "0.0 | \n", "
| \n", " | x | \n", "p | \n", "t | \n", "c | \n", "
|---|---|---|---|---|
| 604 | \n", "[0, 50, 0] | \n", "[2.5, 0.0, 0.0] | \n", "0.00 | \n", "0.5 | \n", "
| 605 | \n", "[0.0495, 49.99902, 0.0] | \n", "[2.475, -0.049, 0.0] | \n", "0.01 | \n", "0.5 | \n", "
| 606 | \n", "[0.09850500000000001, 49.9970698, 0.0] | \n", "[2.45025, -0.09751, 0.0] | \n", "0.02 | \n", "0.5 | \n", "
| \n", " | x | \n", "p | \n", "t | \n", "c | \n", "
|---|---|---|---|---|
| 0 | \n", "50.000000 | \n", "0.000000 | \n", "0.00 | \n", "0.1 | \n", "
| 1 | \n", "49.999020 | \n", "-0.049000 | \n", "0.01 | \n", "0.1 | \n", "
| 2 | \n", "49.997058 | \n", "-0.098098 | \n", "0.02 | \n", "0.1 | \n", "
| 3 | \n", "49.994112 | \n", "-0.147294 | \n", "0.03 | \n", "0.1 | \n", "
| 4 | \n", "49.990180 | \n", "-0.196589 | \n", "0.04 | \n", "0.1 | \n", "
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 284 | \n", "1.671112 | \n", "-18.711466 | \n", "2.84 | \n", "0.1 | \n", "
| 285 | \n", "1.295154 | \n", "-18.797889 | \n", "2.85 | \n", "0.1 | \n", "
| 286 | \n", "0.917464 | \n", "-18.884485 | \n", "2.86 | \n", "0.1 | \n", "
| 287 | \n", "0.538039 | \n", "-18.971254 | \n", "2.87 | \n", "0.1 | \n", "
| 288 | \n", "0.156875 | \n", "-19.058196 | \n", "2.88 | \n", "0.1 | \n", "
289 rows × 4 columns
\n", "