7. Differential Equations

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.

See:

• http://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture21.pdf

• http://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture22.pdf

• Infectious Disease Modelling:

  • Understanding the models that are used to model Coronavirus: Explaining the background and deriving the formulas of the SIR model from scratch. Coding and visualizing the model in Python.

  • Beyond the Basic SIR Model

---

• First Order Methods

  • Euler’s method

  • Example 1

• High Order Methods

  • Second-order Runge-Kutta methods

  • Example 2

  • Fourth-order Runge-Kutta method

• Two-Point Boundary Value Problems

  • Example 3

• SIR Model

---

%pylab inline

Populating the interactive namespace from numpy and matplotlib

import numpy as np
## JSAnimation import available at https://github.com/jakevdp/JSAnimation
#from JSAnimation import IPython_display
from matplotlib import animation

---

7.1. Physical motivation

7.1.1. Momentum

In absence a forces a body moves freely to a constant velocity (See figure)

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

(7.1)\[\begin{align} \boldsymbol{p}=&\gamma m \mathbf{v}\,, &\text{where: } \gamma=&\frac{1}{\sqrt{1-{|\boldsymbol{v}|^2}/{c^2}}}\,, \end{align}\]

and \(c\approx 3\times 10^8\ \)m/s is the speed of light.

If \(\gamma\approx1\), or equivalent, sif \(|\boldsymbol{v}|\ll c\):

(7.2)\[\begin{align} \boldsymbol{p}\approx &m \boldsymbol{v} \end{align}\]

7.1.2. Equation of motion

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

The change of the momentum of particle is equal to the net force acting upon the system times the duration of the interaction

For one instaneous duration \(\Delta t\), this law can be written as $\( \Delta \boldsymbol{p} = \boldsymbol{F}\Delta t,\qquad \Delta t \to 0 \)$ or as the equation of motion, which is the differential equation

(7.3)\[\begin{align} \boldsymbol{F}=&\lim_{t\to 0} \frac{\Delta \boldsymbol{p}}{\Delta t} \\ \boldsymbol{F}=&\frac{\operatorname{d}\boldsymbol{p}}{\operatorname{d} t} \,. \end{align}\]

This equation of motion is of general validity and can be applied numerically to solve any problem.

7.1.3. Constant mass

In the special case of constant mass

(7.4)\[\begin{align} \boldsymbol{F}=&m \boldsymbol{a} \end{align}\]

7.1.4. Example: Drag force

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] $\( F_D=\frac{1}{2}{\rho A C_d} v^2(t)=\frac{1}{2m^2}{\rho A C_d} p^2(t)\,, \)$ where

• \(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

• \(\rho\) is the density of the air

• \(C_d\): Drag coefficient. \(C_d=0.346\) for a baseball

• \(v(t)\): speed of the baseball

• \(p(t)\): momentum of baseball

The equation of motion is

(7.5)\[\begin{align} W-F_D=& \frac{d p}{dt}\\ \frac{d p}{dt}=&m g - \frac{\rho A C_d}{2m^2} p^2(t). \end{align}\]

We see then that, in general $\( \frac{d p}{dt}=f(t,p)\,. \)$

---

7.2. Mathematical background

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.

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.

7.3. Euler’s method

7.3.1. First order differential equations

This method is the most basic of all methods, however, it is useful to understand basic concepts and definitions of differential equations problems.

Suppose we have a well-posed initial-value problem given by:

\[ \frac{dy}{dt}=f(t,y),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \]

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

\[ t_i = a+ i \Delta t,\ \ \ \ \mbox{with}\ \ i=1,2,3,\cdots,N \ \ \mbox{and}\ \Delta t = \frac{b-a}{N}. \]

Here, \(h\) is called the step size of the mesh points.

Now, using the first-order approximation of the derivative that we saw in Numerical Differentiation, we obtain

\[ \frac{dy}{dt}\approx \frac{y(t+\Delta t)-y(t)}{\Delta t} \]

or

\[ \left.\frac{dy}{dt}\right|_{t=t_i}\approx \frac{y(t_{i+1})-y(t_i)}{\Delta t} \]

The original problem can be rewritten as

\[ y_{i+1} = y_i + f(t_i, y_i)\Delta t \pm \frac{\Delta t^2}{2}y^{''}(\xi_i) \]

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''\)

This equation can be solved recursively as we know the initial value \(y_0 = y(a) = \alpha\).

7.3.1.1. Second order differential equations

The formalism of the Euler’s method can be applied for any system of the form:

\[ \frac{dy}{dt}=f(t,y),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \]

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:

\[ \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\]

For this system we have to provide both, the initial value \(y(a)\) and the initial derivative \(y'(a)\).

Now, let’s define a new variable \(w(t) = y'(t)\), the previous problem can be then written as

(7.6)\[\begin{align} \frac{dy}{dt}=&w(t) \\ \frac{dw}{dt}=&f(t,y)-g(t,y)w\, ,\ \ \ a\leq t\leq b\,, \ \ \ \ y(a) = \alpha\ \ \mbox{and}\ w(a) = \beta \end{align}\]

Each of this problem has now the form required by the Euler’s algorithm, and the solution can be computed as:

(7.7)\[\begin{align} w_{i+1}=& w_{i} + [f(t_i,y_i)-g(t_i,y_i)w_i]\Delta t \\ y_{i+1}=&y_i + w_i \Delta t\, ,\ \ \ a\leq t\leq b\,, \ \ \ \ y(a) = \alpha\ \ \mbox{and}\ w(a) = \beta \end{align}\]

7.3.1.2. Euler method for second order differential equations

We can define the column vectors

(7.8)\[\begin{align} \boldsymbol{U}=\begin{bmatrix} U_1\\ U_2\\ \end{bmatrix}=&\begin{bmatrix} y(t)\\ w(t)\\ \end{bmatrix}\,,& \boldsymbol{V}(\boldsymbol{U},t)=\begin{bmatrix} w(t)\\ f(t,y)-g(t,y)w(t)\\ \end{bmatrix}=&\begin{bmatrix} U_2\\ f(t,U_1)-g(t,U_1)U_2\\ \end{bmatrix} \end{align}\]

such that

\[ \frac{dU}{dt}=V(t,U) \]

7.3.2. Newton’s second law of motion

In fact, the equation of motion can be seen as the system of equations

(7.9)\[\begin{align} \frac{\boldsymbol{p}}{m}=&\frac{\operatorname{d}\boldsymbol{x}}{\operatorname{d}t}\\ \boldsymbol{p}=& \gamma m \boldsymbol{v}\\ \boldsymbol{F}=&\frac{\operatorname{d}\boldsymbol{p}}{\operatorname{d} t} \end{align}\]

7.3.3. Example:

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

• Neglecting the air friction

import numpy as np
Δt=0.01 #s
np.arange( 0,3+Δt,Δt)[:2]

array([0.  , 0.01])

(7.10)\[\begin{align} \frac{\boldsymbol{p}}{m}=&\frac{\operatorname{d}\boldsymbol{x}}{\operatorname{d}t}\\ \boldsymbol{p}=& \gamma m \boldsymbol{v}\\ \boldsymbol{F}=&\frac{\operatorname{d}\boldsymbol{p}}{\operatorname{d} t} \end{align}\]

import numpy as np
import pandas as pd
df=pd.DataFrame()
m=.5 ## Kg
g=9.8 ## m/s^2
## intial conditions
x=np.array([0,50,0]) #m
v=np.array([5,0,0]) ## m/s
p=m*v
Δt=0.01 #s
## Analysis for the first 3 s
df=df.append({'x':x,'p':p},ignore_index=True)
Fg=np.array([0,-m*g,0]) #N (Weight)

for t in np.arange( 0,3+Δt,Δt): #s
    p=p+Fg*Δt
    x=x+(p/m)*Δt
    df=df.append({'x':x,'p':p},ignore_index=True)

df[:3]

	x	p
0	[0, 50, 0]	[2.5, 0.0, 0.0]
1	[0.05, 49.99902, 0.0]	[2.5, -0.049, 0.0]
2	[0.1, 49.997060000000005, 0.0]	[2.5, -0.098, 0.0]

df.to_json('mvto.json')

#pd.read_json('mvto.json').x.str[1]

Remember that in Python an string is also a list:

s='abc'

len(s)

3

s[0]

'a'

s[-1]

'c'

We can use the str attribute

#Second component of vector x:
df['x'].str[1][:3]

0    50.00000
1    49.99902
2    49.99706
Name: x, dtype: float64

#py=df.p.str[1]

##               x       y
plt.plot(df.x.str[0],df.x.str[1])
plt.xlabel('$x$ [m]',size=15)
plt.ylabel('$y$ [m]',size=15)
plt.grid()

Consider now the case the movement inside a fluid with a dragging force proportinal to the velocity (low velocity case) $\( \boldsymbol{F}_D=-c \boldsymbol{v} = -c \frac{\boldsymbol{p}}{m}\,, \)\( where \)c$ can take the values 0 (not dragging force), 0.1 y 0.5

import numpy as np
import pandas as pd
df=pd.DataFrame()
m=.5 #kg
g=9.8 #m/s^2
for c in [0,0.1,0.5]:
    x=np.array([0,50,0]) #m
    v=np.array([5,0,0]) #m/s
    p=m*v
    Δt=0.01 #s
    t=0 #s
    df=df.append({'x':x,'p':p,'t':t,'c':c},ignore_index=True)
    
    for t in np.arange( 0,3+Δt,Δt):
        F=np.array([0,-m*g,0])-c*(p/m)
        p=p+F*Δt
        x=x+(p/m)*Δt
        t=t+Δt
        df=df.append({'x':x,'p':p,'t':t,'c':c},ignore_index=True)

df.c.unique()

array([0. , 0.1, 0.5])

Apply a mask upon the DataFrame:

• Example

• Filter c==0

df[df.c==0][:3]

	x	p	t	c
0	[0, 50, 0]	[2.5, 0.0, 0.0]	0.00	0.0
1	[0.05, 49.99902, 0.0]	[2.5, -0.049, 0.0]	0.01	0.0
2	[0.1, 49.997060000000005, 0.0]	[2.5, -0.098, 0.0]	0.02	0.0

df[df.c==0.5][:3]

	x	p	t	c
604	[0, 50, 0]	[2.5, 0.0, 0.0]	0.00	0.5
605	[0.0495, 49.99902, 0.0]	[2.475, -0.049, 0.0]	0.01	0.5
606	[0.09850500000000001, 49.9970698, 0.0]	[2.45025, -0.09751, 0.0]	0.02	0.5

%pylab inline

Populating the interactive namespace from numpy and matplotlib

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.gca(projection='3d')
c=0
ax.plot(df[df.c==c].x.str[0].values, df[df.c==c].x.str[1].values, df[df.c==c].x.str[2].values)
c=0.1
ax.plot(df[df.c==c].x.str[0].values, df[df.c==c].x.str[1].values, df[df.c==c].x.str[2].values)
c=0.5
ax.plot(df[df.c==c].x.str[0].values, df[df.c==c].x.str[1].values, df[df.c==c].x.str[2].values)

ax.view_init(80, 250)
plt.xlabel('$x$ [m]',size=10)
plt.ylabel('$y$ [m]',size=10)

/tmp/ipykernel_61048/2642031594.py:2: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')

Text(0.5, 0, '$y$ [m]')

c=0
plt.plot(df[df.c==c].x.str[0], df[df.c==c].x.str[1],label=r'$c=%s$ ' %c)
c=0.1
plt.plot(df[df.c==c].x.str[0], df[df.c==c].x.str[1],label=r'$c=%s$ ' %c)
c=0.5
plt.plot(df[df.c==c].x.str[0], df[df.c==c].x.str[1],label=r'$c=%s$ ' %c)
plt.xlabel('$x$ [m]',size=15)
plt.ylabel('$y$ [m]',size=15)
plt.legend(loc='best')
plt.grid()

Activity: Find the time to reach the position y=0 in each case. By using an algorithm to find roots

Activity: Repeat the previous analysis for a baseball in the air with a dragging force proportinal to the square of the velocity. The diameter of the baseball is 7.5cm

7.3.3.1. Example. Spring

In order to apply this, let’s assume a simple mass-spring.

The equation of motion according to Newton’s second law is

\[ F = -kx \]

Using the previous results, we can rewrite this as:

\[ \frac{d p}{dt} = -k x \]

\[ \frac{dx}{dt} = \frac{p}{m}\]

And the equivalent Euler system is

(7.11)\[\begin{align} v_{i+1}=& v_{i} - \Delta t\frac{k}{m}x_i \\ x_{i+1}=&x_i + v_i \Delta t\, ,\ \ \ x(0) = x_0\ \ \mbox{and}\ v(0) = v_0 \end{align}\]

7.3.4. Activity

**1.** Using the initial conditions $x(0) = 0$ and $v(0) = 3$, solve the previous system. Plot the solutions $x(t)$ and $y(t)$ and compare with real solutions. Furthermore, calculate the total energy of the system. What can you conclude about the behaviour of the energy? Does it make any sense? **2.** Using the same reasoning, derive the equations for a simple pendulum. Compare the solution for small oscillations with the analytic one. What happens when increasing the initial amplitude of the movement?

---

(7.12)\[\begin{align} \frac{\boldsymbol{p}}{m}=&\frac{\operatorname{d}\boldsymbol{x}}{\operatorname{d}t}\\ \boldsymbol{p}=& \gamma m \boldsymbol{v}\\ \boldsymbol{F}=&\frac{\operatorname{d}\boldsymbol{p}}{\operatorname{d} t} \end{align}\]

1D $\( x_{i+1}=x_i+(p/m)\Delta t\)\( \)\( p_{i+1}-p_i= F\Delta t\)\( \)\( p_{i+1}=p_i+F \Delta t \)$

import numpy as np
import pandas as pd
df=pd.DataFrame()
m=.5 #kg
g=9.8 #m/s^2
c=0.1
x=50 #m
v=0 #m/s
p=m*v
Δt=0.01 #s
t=0 #s
df=df.append({'x':x,'p':p,'t':t,'c':c},ignore_index=True)
    
for t in np.arange( 0,3+Δt,Δt):
        F=-m*g+c*(p/m)
        p=p+F*Δt
        x=x+(p/m)*Δt
        t=t+Δt
        df=df.append({'x':x,'p':p,'t':t,'c':c},ignore_index=True)

df[df['x']>=0]

	x	p	t	c
0	50.000000	0.000000	0.00	0.1
1	49.999020	-0.049000	0.01	0.1
2	49.997058	-0.098098	0.02	0.1
3	49.994112	-0.147294	0.03	0.1
4	49.990180	-0.196589	0.04	0.1
...	...	...	...	...
284	1.671112	-18.711466	2.84	0.1
285	1.295154	-18.797889	2.85	0.1
286	0.917464	-18.884485	2.86	0.1
287	0.538039	-18.971254	2.87	0.1
288	0.156875	-19.058196	2.88	0.1

289 rows × 4 columns

7.3.5. Implementation in Scipy

%pylab inline

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

/usr/local/lib/python3.9/dist-packages/IPython/core/magics/pylab.py:162: UserWarning: pylab import has clobbered these variables: ['f']
`%matplotlib` prevents importing * from pylab and numpy
  warn("pylab import has clobbered these variables: %s"  % clobbered +

import scipy.integrate as integrate
import numpy as np

7.3.5.1. First order ordinary differential equations

integrate.odeint: Integrate a system of ordinary differential equations

Uso básico

integrate.odeint(func,α,t)

Solves the initial value problem for stiff or non-stiff systems of first order ordinary differential equations:

\[ \frac{dy}{dt}=f(y,t) \]

where \(y\) can be a vector. \(a<t<b\) and \(y(a)=\alpha\)

7.3.5.2. Example

Consider the following differential equation

\[\frac{dy(t)}{dt}=-k y(t),\]

where \(k\) is constant and with initial condition \(y(0)=y_0=5\)

def f(y,t): #t is mandatory even if no explicit 
    k = 0.3
    return -k * y

For a first evolution after one step of \(\Delta t=0.1\): with initial condition \(y(0)=y_0=5\)

## initial condition
y0 = 5
t=[0,0.1]
integrate.odeint(f,y0,t)

array([[5.        ],
       [4.85222774]])

while the evolution from 0 to 20 s is

y=y0+f(y0,t)*0.1
y

4.85

## time points
t = np.linspace(0,20)

## solve ODE
y = integrate.odeint(f,y0,t)

## plot results
plt.plot(t,y)
plt.xlabel('$t$',size=15)
plt.ylabel('$y(t)$',size=15)
plt.grid()
plt.show()

Compare with the analytical solution $\(y(t)=y_0\operatorname{e}^{-kt}\)$

Since $\( \frac{dy}{dt}=\left[\frac{d (-kt)}{dt}\right] y_0\operatorname{e}^{-kt}=-k \left(y_0\operatorname{e}^{-kt}\right)=-k\,y(t) \)$

## time points
t = np.linspace(0,20)

#k = 0.3 solve ODE
y0=5
y = integrate.odeint(f,y0,t)

k=0.3
## plot results
plt.plot(t,y,'c-',label='ODE solution')
plt.plot(t,y0*np.exp(-k*t),'k:',label='Analytical solution')
plt.xlabel('$t$',size=15)
plt.ylabel('$y(t)$',size=15)
plt.grid()
plt.legend(loc='best')
plt.show()

7.3.5.3. Second order differential equations

Although first-order schemes like Euler’s method are illustrative and allow a good understanding of the numerical problem, real applications cannot be dealt with them, instead more precise and accurate high-order methods must be invoked. In this section we shall cover a well-known family of numerical integrators, the Runge-Kutta methods.

import scipy.integrate as integrate

7.3.5.4. Some examples

1. http://www.southampton.ac.uk/~fangohr/teaching/python/book/Python-for-Computational-Science-and-Engineering.pdf

   • http://www.southampton.ac.uk/~fangohr/teaching/python/book/ipynb/

2. http://sam-dolan.staff.shef.ac.uk/mas212/

   • http://sam-dolan.staff.shef.ac.uk/mas212/notebooks/ODE_Example.ipynb

3. https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

4. http://csc.ucdavis.edu/~cmg/Group/readings/pythonissue_3of4.pdf

5. http://www.tau.ac.il/~kineret/amit/scipy_tutorial/#x1-210004.4

6. https://www.udacity.com/course/differential-equations-in-action–cs222

7. https://github.com/robclewley/pydstool

   • http://www2.gsu.edu/~matrhc/PyDSTool.htm

8. https://github.com/JuliaDiffEq

Example: From http://sam-dolan.staff.shef.ac.uk/mas212/

As explained before, we need to write a second order ordinary differential equations in terms of first order matricial ordinary differential equation. In terms of a parameter, \(x\), this implay to have a column vector $\( U=\begin{bmatrix} y\\ z\\ \end{bmatrix} \)\(4 \)\( \frac{dU}{dt}=V(U,x). \)$

Suppose we have a second-order ODE such as a damped simple harmonic motion equation, with parameter \(x\) (instead \(t\)) $\( \quad y'' + 2 y' + 2 y = \cos(2x), \quad \quad y(0) = 0, \; y'(0) = 0 \)\( We can turn this into two first-order equations by defining a new depedent variable. For example, \)\( \quad z \equiv y' \quad \Rightarrow \quad z' + 2 z + 2y = \cos(2x), \quad z(0)=y(0) = 0. \)\( So that \)\( z' =- 2 z - 2y + \cos(2x), \quad z(0)=y(0) = 0. \)\( where \)\( \frac{dU}{dx}=\begin{bmatrix}y'\\z'\end{bmatrix} \)$ We can solve this system of ODEs using “odeint” with lists, as follows:

Let $\( U=\begin{bmatrix} U_0\\ U_1 \end{bmatrix}=\begin{bmatrix} y\\ z \end{bmatrix}. \)$ Therefore

(7.13)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} x} \begin{bmatrix} y\\ z\\ \end{bmatrix}=& \begin{bmatrix} z\\ -2z-2y+\cos(2x)\\ \end{bmatrix}\\ \frac{\operatorname{d}}{\operatorname{d} x} U=& \begin{bmatrix} U_1\\ -2U_1-2U_0+\cos(2x)\\ \end{bmatrix} \end{align}\]

Implementation by using only \(U\)

def dU_dx(U, x):
    '''
    Here U is a vector such that y=U[0] and z=U[1]. 
    This function should return [y', z']
    '''    
    return [      U[1], 
            -2*U[1] - 2*U[0] + np.cos(2*x)]

U0=[0,0] ## → x=0

x→0.1

integrate.odeint(dU_dx, U0, [0,0.1])

array([[0.        , 0.        ],
       [0.00465902, 0.08970031]])

array(
        [[0.        , 0.         ],#→ U(0)
         [0.00465902, 0.08970031]] #→ U(0.1)=[y(0.1),z(0.1)=y'(0.1)]
     )

array([[0.        , 0.        ],
       [0.00465902, 0.08970031]])

U0 = [0, 
      0]
xs = np.linspace(0, 50, 200)
Us = integrate.odeint(dU_dx, U0, xs)

Us[:3]

array([[0.        , 0.        ],
       [0.02601336, 0.1841925 ],
       [0.08083294, 0.229454  ]])

Us.shape

(200, 2)

The first column, Us[:,0] → y

Us[:,0][:3]

array([0.        , 0.02601336, 0.08083294])

ys = Us[:,0] #First column is extracted

plt.xlabel("$x$",size=15)
plt.ylabel("$y$",size=15)
plt.title("Damped harmonic oscillator")
plt.plot(xs,ys);

zs = Us[:,1]

plt.xlabel("$x$",size=15)
plt.ylabel("$z=y'$",size=15)
plt.title("Damped harmonic oscillator")
plt.plot(xs,zs);

Implementation by using explicitly \(y(x)\) and \(z(x)\)

(7.14)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} x} \begin{bmatrix} y\\ z\\ \end{bmatrix}=& \begin{bmatrix} z\\ -2z-2y+\cos(2x)\\ \end{bmatrix} \end{align}\]

def dU_dx(U, x):
    '''
    Here U is a vector such that y=U[0] and z=U[1]. 
    This function should return [y', z']
    '''    
    y,z=U ## y → U[0]; z → U[1] 
    return [        z, 
            -2*z - 2*y + np.cos(2*x)]

U0 = [0, 0]
xs = np.linspace(0, 50, 200)
Us = integrate.odeint(dU_dx, U0, xs)

ys = Us[:,0]

plt.xlabel("$x$",size=15)
plt.ylabel("$y$",size=15)
plt.title("Damped harmonic oscillator")
plt.plot(xs,ys);

Activity: Apply the previous example to the problem of parabolic motion with air friction:

(7.15)\[\begin{align} \frac{\boldsymbol{p}}{m}=&\frac{\operatorname{d}\boldsymbol{x}}{\operatorname{d}t}\\ \boldsymbol{F}=&\frac{\operatorname{d}\boldsymbol{p}}{\operatorname{d} t} \end{align}\]

(7.16)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} t} \begin{bmatrix} \boldsymbol{x}\\ \boldsymbol{p}\\ \end{bmatrix}=& \begin{bmatrix} \boldsymbol{p}/m\\ -m\boldsymbol{g}-c \boldsymbol{p}/m \\ \end{bmatrix}\,, \end{align}\]

where \(\boldsymbol{g}=(0,g)\), and such that in two dimensions

\[\begin{split} \boldsymbol{U}=\begin{bmatrix} {U}_0\\ {U}_1\\ {U}_2\\ {U}_3\\ \end{bmatrix}=\begin{bmatrix} {x}\\ y\\ p_x\\ p_y\\ \end{bmatrix} \end{split}\]

(7.17)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} t} \boldsymbol{U}=& \begin{bmatrix} {U}_2/m\\ {U}_3/m\\ -c{U}_2/m\\ -m{g}-c{U}_3/m\\ \end{bmatrix} \end{align}\]

def dU_dt(U, t,c=0.,m=.5,g=9.8):
    '''
    Here U is a vector such that y=U[0] and z=U[1]. 
    This function should return [y', z']
    '''
    return [U[2]/m,
            U[3]/m,
            -c*U[2]/m, 
            -m*g-c*U[3]/m]

m=0.5
## p_x(0)=5*m → v(0)=5 m/s
U0 = [0, 50,m*5,0]
t = np.linspace(0, 3, 200)

isinstance((0.1),float)

True

isinstance((0.1,),tuple)

True

Us = integrate.odeint(dU_dt, U0, t )  
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=r'$c=0$')
Us = integrate.odeint(dU_dt, U0, t,args=(0.1,) )  
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=r'$c=0.1$')
Us = integrate.odeint(dU_dt, U0, t,args=(0.5,) )  
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=r'$c=0,5$')
plt.xlabel('$x$ [m]',size=15)
plt.ylabel('$y$ [m]',size=15)
plt.legend(loc='best')
plt.grid()

Activity: Implemet directly:

(7.18)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} t} \begin{bmatrix} \boldsymbol{x}\\ \boldsymbol{p}\\ \end{bmatrix}=& \begin{bmatrix} \boldsymbol{p}/m\\ -m\boldsymbol{g}-c \boldsymbol{p}/m \\ \end{bmatrix}\,, \end{align}\]

as

(7.19)\[\begin{align} \frac{\operatorname{d}}{\operatorname{d} t} \begin{bmatrix} x\\ y\\ p_x\\ p_y \end{bmatrix}=& \begin{bmatrix} p_x/m\\ p_y/m\\ -c p_x/m \\ -mg-c p_y/m \\ \end{bmatrix}\,, \end{align}\]

def dU_dt(U, t,c=0.,m=.5,g=9.8):
    '''
    ...
    '''
    x,y,px,py=U
    return [px/m,
            py/m,
            -c*px/m, 
            -m*g-c*py/m]

m=0.5
U0 = [0, 50,5*m,0]
t = np.linspace(0, 3, 200)

c=0
Us = integrate.odeint(dU_dt, U0, t,args=(c,) ) ## c = 0.5 
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=f'c={c}')


c=0.1
Us = integrate.odeint(dU_dt, U0, t,args=(c,) ) ## c = 0.5 
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=f'c={c}')

c=0.5
Us = integrate.odeint(dU_dt, U0, t,args=(c,) ) ## c = 0.5 
xs = Us[:,0]
ys = Us[:,1]
plt.plot(xs,ys,label=f'c={c}')

plt.xlabel('$x$ [m]',size=15 )
plt.ylabel('$y$ [m]',size=15 )
plt.legend(loc='best')
plt.grid()

%pylab inline

Populating the interactive namespace from numpy and matplotlib

7.3.6. Example 2

In this example we are going to use the Scipy implementation for mapping the phase space of a pendulum.

The equations of the pendulum are given by:

\[ \frac{d\theta}{dt} = \omega \]

\[ \frac{d\omega}{dt} = -\frac{g}{l}\sin \theta \]

(7.20)\[\begin{align} \frac{d}{dt} U=\begin{pmatrix}\omega\\ -\dfrac{g}{l}\sin\theta \end{pmatrix}, \end{align}\]

where

(7.21)\[\begin{align} U=\begin{pmatrix}\theta\\ \omega \end{pmatrix}. \end{align}\]

from scipy import integrate
def dU_dt(U, t,l=1,g=9.8):
    '''
    Here U is a vector such that   θ=U[0] and ω=U[1]. 
    This function should return [θ', ω']
    '''    
    θ,ω=U
    return [ω, -g/l*np.sin( θ ) ]

tmax = 6*np.pi #s
omega_max = 8 #rad/s

Nic = 1000
#Maxim angular velocity

theta0s = np.random.uniform(-4*np.pi,4*np.pi,Nic)
omega0s = np.random.uniform(-omega_max,omega_max,Nic)

TAREA: Generar un rango aleatorio uniforme de varios ordenes de magnitud. En particular, 1000 números aleatorios entre \(3\times 10^{-6}\) hasta \(5\times 10^{4}\)

theta0s[0]

1.3936715764887904

U0=[theta0s[0],omega0s[0]]
U0

[1.3936715764887904, 7.007902033513252]

t=np.linspace(0,tmax,400)
Us=integrate.odeint(dU_dt,U0,t)

plt.plot(Us[:,0],Us[:,1],lw = 1, color = "black" )
plt.xlim(-10,10)
plt.ylim(-10,10)
plt.xlabel(r'$\theta$ [rad]',size=15)
plt.ylabel(r'$\omega$ [rad/s]',size=15)

Text(0, 0.5, '$\\omega$ [rad/s]')

list( zip([1,2],[3,4]) )

[(1, 3), (2, 4)]

j=0
plt.figure( figsize = (8,5) )
for theta0, omega0 in zip(theta0s, omega0s):
    t=np.linspace(0,tmax,400)
    U0=[theta0,omega0]
    Us=integrate.odeint(dU_dt,U0,t)
    plt.plot(Us[:,0]/np.pi,Us[:,1],lw = 0.1, color = "black" )
    
#Format of figure
plt.xlabel( "$\\theta/\pi$", fontsize = 18 )
plt.ylabel( "$\omega$ [rad/s]", fontsize = 18 )
plt.xlim( (-3, 3) )
plt.ylim( (-omega_max, omega_max) )
plt.title( "Phase space of a pendulum" )
plt.grid(1)    

The nearly closed curves around (0,0) represent the regular small swings of the pendulum near its rest position. The oscillatory curves up and down of the closed curves represent the movement when the pendulum start at \(\theta=0\) but with high enough angular speed such that the pendulum goes all the way around. Of course, its angular speed will slow down on the way up but then it will speed up on the way down again. In the absence of friction, it just keeps spinning around indefinitely. The counterclockwise motions of the pendulum of this kind are shown in the graph by the wavy lines at the top with positive angular speed, while the curves on the bottom which go from right to left represent clockwise rotations. The phase space have a periodicity of \(2\pi\).

Small oscillation

t=np.linspace(0,tmax,400)
Us=integrate.odeint(dU_dt,[0,2],t)
plt.plot(Us[:,0],Us[:,1],lw = 1, color = "black" )
plt.xlim(-8,8)
plt.ylim(-10,10)

(-10.0, 10.0)

All around

t=np.linspace(0,tmax,400)
Us=integrate.odeint(dU_dt,[0,7],t)
plt.plot(Us[:,0],Us[:,1],lw = 1, color = "black" )
plt.xlim(-8,8)
plt.ylim(-10,10)

(-10.0, 10.0)

t=np.linspace(0,tmax,400)
Us=integrate.odeint(dU_dt,[0,-7],t)
plt.plot(Us[:,0],Us[:,1],lw = 1, color = "black" )
plt.xlim(-8,8)
plt.ylim(-10,10)

(-10.0, 10.0)

plt.plot(Us[:,0],Us[:,1],lw = 1, color = "black" )
plt.xlim(-8,8)
plt.ylim(-10,10)

(-10.0, 10.0)

Activity: Check the anlytical expression for the period of a pendulum

---

7.4. Appendix

7.4.1. Activity

Using the previous example, explore the phase space of a simple oscillator and a damped pendulum.

7.4.2. Fourth-order Runge-Kutta method

Finally, the most used general purpose method is the fourth-order Runge-Kutta scheme. Its derivation follows the same previous reasoning, however the procedure is rather long and it makes no sense to reprouce it here. Instead, we will give the direct algorithm:

Let’s assume again a problem of the form:

\[ \frac{dy}{dt}=y'=f(t,y),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \]

The Runge-Kutta-4 (RK4) method allows us to predict the solution at the time \(t+h\) as:

\[ y(t+h) = y(t) + \frac{1}{6}( \mathbf{K}_0 + 2\mathbf{K}_1 + 2\mathbf{K}_2 + \mathbf{K}_3 ) \]

where:

\[ \mathbf{K}_0 = hf(t,y)\]

\[ \mathbf{K}_1 = hf\left( t + \frac{h}{2},y + \frac{\mathbf{K}_0}{2}\right)\]

\[ \mathbf{K}_2 = hf\left( t + \frac{h}{2},y + \frac{\mathbf{K}_1}{2}\right)\]

\[ \mathbf{K}_3 = hf\left( t + h,y + \mathbf{K}_2\right)\]

7.4.3. Activity

The Lorenz attractor is a common set of differential equations of some models of terrestrial atmosphere studies. It is historically known as one of the first system to exhibit deterministic caos. The equations are:

\[ \frac{dx}{dt} = a(y-x) \]

\[ \frac{dy}{dt} = x(b-z)-y \]

\[ \frac{dz}{dt} = xy-cz \]

with \(a = 10\), \(b=28\) and \(c = 8/3\) the solution shows periodic orbits.

Write a routine that gives a step of RK4 and integrate the previous system. Plot the resulting 3D solution $(x,y,z)$.

7.4.4. Second-order Runge-Kutta methods

For this method, let’s assume a problem of the form:

\[ \frac{dy}{dt}=y'=f(t,y),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \]

Now, we want to know the solution in the next timestep, i.e. \(y(t+h)\). For this, we propose the next solution:

\[ y(t+h) = y(t) + c_0 f(t,y)h + c_1f[ t+ph, y+qhf(t,y) ]h \]

determining the coefficients \(c_0, c_1, q\) and \(p\), we will have the complete approximated solution of the problem.

One way to determine them is by comparing with the taylor expansion around \(t\)

\[ y(t+h) = y(t) + f(t,y)h + \frac{1}{2}\left( \frac{\partial f}{\partial t} + \frac{\partial f}{\partial y} \right)h^2 + \mathcal{O}(h^3) \]

Now, we can also expand the function \(f[ t+ph, y+qhf(t,y) ]\) around the point \((t,y)\), yielding:

\[ f[ t+ph, y+qhf(t,y) ] = f(t,y) + \frac{\partial f}{\partial t}ph + \frac{\partial f}{\partial y}qh + \mathcal{O}(h^2) \]

Replacing this into the original expression:

\[ y(t+h) = y(t) + c_0 f(t,y)h + c_1\left[ f(t,y) + \frac{\partial f}{\partial t}ph + \frac{\partial f}{\partial y}qh \right]h + \mathcal{O}(h^3)\]

ordering the terms we obtain:

\[ y(t+h) = y(t) + (c_0+c_1) f(t,y)h + c_1\left[ \frac{\partial f}{\partial t}p + \frac{\partial f}{\partial y}q \right]h^2 + \mathcal{O}(h^3)\]

Equalling the next conditions are obtained:

\[ c_0 + c_1 =1 \ \ \ c_1p=\frac{1}{2}\ \ \ c_1q = \frac{1}{2} \]

This set of equations are undetermined so there are several solutions, each one yielding a different method:

\[\begin{split} \matrix{ c_0 = 0 & c_1=1 & p = 1/2 & q = 1/2 & \mbox{Modified Euler's Method} \\ c_0 = 1/2 & c_1=1/2 & p = 1 & q = 1 & \mbox{Heun's Method} \\ c_0 = 1/3 & c_1=2/3 & p = 3/4 & q = 3/4 & \mbox{Ralston's Method} } \end{split}\]

The algorithm is then:

\[ y(t+h) = y(t) + (c_1+c_2)\mathbf{K}_1 \]

with

\[ \mathbf{K}_1 = hf( t+pt, y+q\mathbf{K}_0 ) \]

\[ \mathbf{K}_0 = hf(t,y) \]

All these methods constitute the second-order Runge-Kutta methods.

7.5. Two-Point Boundary Value Problems

Up to now we have solved initial value problems of the form:

\[ \frac{dy}{dt}=y'=f(t,y),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \]

Second order equations can be similarly planted as

\[ \frac{d^2y}{dt^2}=y''=f(t,y,y'),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \ \ \ y'(a) = u \]

This type of systems can be readily solved by defining the auxiliar variable \(w = y'\), turning it into a first order system of equations.

Now, we shall solve two-point boundary problem, where we have two conditions on the solution \(y(t)\) instead of having the function and its derivative at some initial point, i.e.

\[ \frac{d^2y}{dt^2}=y''=f(t,y,y'),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \ \ \ y(b) = \beta \]

In spite of its similar form to the initial value problem, two-point boundary problems pose a increased difficulty for numerical methods. The main reason of this is the iterative procedure performed by numerical approaches, where from an initial point, further points are found. Trying to fit two different values at different points implies then a continuous readjusting of the solution.

A common way to solve these problems is by turning them into a initial-value problem

\[ \frac{d^2y}{dt^2}=y''=f(t,y,y'),\ \ \ a\leq t\leq b, \ \ \ \ y(a) = \alpha \ \ \ y'(a) = u \]

Let’s suppose some choice of \(u\), integrating by using some of the previous methods, we obtain the final boundary condition \(y(b)=\theta\). If the produced value is not the one we wanted from our initial problem, we try another value \(u\). This can be repeated until we get a reasonable approach to \(y(b)=\beta\). This method works fine, but it is so expensive and terribly inefficient.

Note when we change \(u\), the final boundary value also change, so we can assume \(y(b) = \theta\). The solution to the problem can be thought then as a root-finding problem:

\[ y(b) = \theta(u) = \beta \]

or

\( r(u) \equiv \theta(u) - \beta = 0 \)

where \(r(u)\) is the residual function. This problem can be thus solved using some of the methods previously seen for the root-finding problem.

7.5.1. Example 3

A very simplified model of interior solid planets consists of a set of spherically symmetric radial layers, where the next properties are computed: density \(\rho(r)\), enclosed mass \(m(r)\), gravity \(g(r)\) and pressure \(P(r)\). Each of these properties are assumed to be only radially dependent. The set of equations that rules the planetary interior is:

Hydrostatic equilibrium equation

\[\frac{dP}{dr} = -\rho(r)g(r)\]

Adams-Williamson equation

\[\frac{dg}{dr} = 4\pi G\rho(r) - \frac{2Gm(r)}{r^3}\]

Continuity equation

\[\frac{dm}{dr} = 4\pi r^2 \rho(r)\]

Equation of state

\[\frac{d\rho}{dr} = -\frac{\rho(r)^2g(r)}{K_s}\]

For accurate results the term \(K_s\), called the adiabatic bulk modulus, is temperature and radii dependent. However, for the sake of simplicity we shall assume a constant value.

Solving simultaneously the previous set of equations, we can find the complete internal structure of a planet.

We have four functions to be determined and four equations, so the problem is solvable. It is only necessary to provide a set of boundary conditions of the form:

\[ \rho(R) = \rho_{surf},\ \ \ m(R) = M_p, \ \ \ g(R) = g_{surf},\ \ \ P(R) = P_{atm} \]

where \(R\) is the planetary radius, \(\rho_{surf}\) the surface density, \(M_p\) the mass of the planet, \(g_{surf}\) the surface gravity and \(P_{atm}\) the atmospheric pressure. However, there is a problem, we do not know the planetary radius \(R\), so an extra condition is required to determine this value. This is reached through the physical condition \(m(0) = 0\), this is, the enclosed mass at a radius \(r = 0\) (center of the planet) must be 0.

The two-value boundary nature of this problem lies then in fitting the mass function at \(m(R) = M_p\) and at \(m(0) = 0\). To do so, let’s call the residual mass \(m(0) = M_r\). This value should depend on the chosen radius \(R\), so a large value would imply a mass defect \(M_r(R)<0\), and a small value a mass excess \(M_r(0)>0\). The problem is then solving the radius \(R\) for which \(m(0) = M_r(R) = 0\). This can be done by using the bisection method.

For this problem, we are going to assume an one-layer planet made of perovskite, so \(K_s \approx 200\ GPa\). A planet mass equal to earth, so \(M_p = 5.97219 \times 10^{24}\ kg\), a surface density \(\rho_{surf} = 3500\ kg/m^3\) and a atmospheric pressure of \(P_{atm} = 1\ atm = 1\times 10^5\ Pa\).

7.5.2. Implementation of the pendulum phase space

#RK2 integrator
def RK2_step( f, y, t, h ):
    #Creating solutions
    K0 = h*f(t, y)
    K1 = h*f(t + 0.5*h, y + 0.5*K0)
    y1 = y + K1
    #Returning solution
    return y1

The phase space of a dynamical system is a space in which all the possible states of that system are univocally represented. For the case of the pendulum, a complete state of the system is given by the set \((\theta, \omega)\), so its phase space is two-dimensional. In order to explore all the possible states, we are going to generate a set of initial conditions and integrate them.

#========================================================
#Defining parameters
#========================================================
#Gravity
g = 9.8
#Pendulum's lenght
l = 1.0

#Number of initial conditions
Nic = 1000
#Maxim angular velocity
omega_max = 8

#Maxim time of integration
tmax = 6*np.pi
#Timestep
h = 0.01

#========================================================
#Dynamical function of the system
#========================================================
def function( t, y ):
    #Using the convention y = [theta, omega]
    theta = y[0]
    omega = y[1]
    #Derivatives
    dtheta = omega
    domega = -g/l*np.sin(theta)
    return np.array([dtheta, domega])

#========================================================
#Generating set of initial conditions
#========================================================
theta0s = -4*np.pi + np.random.random(Nic)*8*np.pi
omega0s = -omega_max + np.random.random(Nic)*2*omega_max

#========================================================
#Integrating and plotting the solution for each IC
#========================================================
#Setting figure
plt.figure( figsize = (8,5) )
jj=0
for theta0, omega0 in zip(theta0s, omega0s):
    #Arrays for solutions
    time = [0,]
    theta = [theta0,]
    omega = [omega0,]
    for i, t in zip(range(int(tmax/h)), np.arange( 0, tmax, h )):
        #Building current condition
        y = [theta[i], omega[i]]
        #Integrating the system
        thetai, omegai = RK2_step( function, y, t, h )
        #Appending new components
        theta.append( thetai )
        omega.append( omegai )
        time.append( t )
        #if i==10:
        ##    break
    #Plotting solution
    plt.plot( theta, omega, lw=0.1, color = "blue" )
    if jj==50:
        break
    jj=jj+1

#Format of figure
plt.xlabel( "$\\theta$", fontsize = 18 )
plt.ylabel( "$\omega$", fontsize = 18 )
plt.xlim( (-2*np.pi-2, 2*np.pi+2) )
plt.ylim( (-omega_max-3, omega_max+3) )
plt.title( "Phase space of a pendulum" )
plt.grid(1)