Marcos Duarte
Laboratory of Biomechanics and Motor Control (http://demotu.org/)
Federal University of ABC, Brazil
An ordinary differential equation (ODE) is an equation containing a function of one independent variable and its derivatives.
Solve an ODE is finding such a function whose derivatives satisfy the equation. The order of an ODE refers to the order of the derivatives; e.g., a first order ODE has only first derivatives. A linear ODE has only linear terms for the function of one independent variable and in general its solution can be obtained analytically. By contrast, a nonlinear ODE doesn't have an exact analytical solution and it has to be solved by numerical methods. The equation is referred as partial differential equation when contains a function of more than one independent variable and its derivatives.
A simple and well known example of ODE is Newton's second law of motion:
$$ m\frac{\mathrm{d}^2 \mathbf{x}}{\mathrm{d}t^2}(t) = \mathbf{F} $$$\mathbf{x}$ is the function with a derivative and $t$ is the independent variable. Note that the force, $\mathbf{F}$, can be constant (e.g., the gravitational force) or a function of position, $\mathbf{F}(\mathbf{x}(t))$, (e.g., the force of a spring) or a function of other quantity. If $\mathbf{F}$ is constant or a linear function of $\mathbf{x}$, this equation is a second-order linear ODE.
A first-order ODE has the general form:
$$ \frac{\mathrm{d} y}{\mathrm{d} x} = f(x, y) $$Where $f(x, y)$ is an expression for the derivative of $y$ that can be evaluated given $x$ and $y$. When $f(x, y)$ is linear w.r.t. $y$, the equation is a first-order linear ODE which can be written in the form:
$$ \frac{\mathrm{d} y}{\mathrm{d} x} + P(x)y = Q(x) $$When an ODE can't be solved analytically, usually because it's nonlinear, numerical methods are used, a procedure also referred as numerical integration (Downey, 2011; Kitchin, 2013; Kiusalaas, 2013; Wikipedia). In numerical methods, a first-order differential equation can be solved as an Initial Value Problem (IVP) of the form:
$$ \dot{y}(t) = f(t, y(t)), \quad y(t_0) = y_0 $$In numerical methods, a higher-order ODE is usually transformed into a system of first-order ODE and then this system is solved using numerical integration.
The most simple method to solve an ODE is using the Euler method.
First, the derivative of $y$ is approximated by:
Where $h$ is the step size.
Rearranging the equation above:
And replacing $\dot{y}(t)$:
$$ y(t+h) \approx y(t) +hf(t, y(t)) $$The ODE then can be solved starting at $t_0$, which has a known value for $y_0$:
$$ y(t+h) \approx y_0 + hf(t_0, y_0) $$And using the equation recursively for a sequence of values for $t$ $(t_0, t_0+h, t_0+2h, ...)$:
$$ y_{n+1} = y_n + hf(t_n, y_n) $$This is the Euler method to solve an ODE with a known initial value.
There are other methods for solving an ODE. One family of methods, usually more accurate, uses more points in the interval $[t_n,t_{n+1}]$ and are known as Runge–Kutta methods. In the Python ecosystem, Runge–Kutta methods are available using the scipy.integrate.ode
library of numeric integrators. The library scipy.integrate.odeint
has other popular integrator known as lsoda
, from the FORTRAN library odepack.
Consider a football ball kicked up from an initial height $y_0$ and with initial velocity $v_0$. Determine the equation of motion of the ball in the vertical direction.
Neglecting the air resistance, Newton's second law of motion applied to this problem for the instants the ball is in the air gives:
$$ m\frac{\mathrm{d}^2 y}{\mathrm{d}t^2} = -mg $$Consider $g=9.8m/s^2$, $y_0(t_0=0)=1m$, and $v_0(t_0=0)=20m/s$.
We know the analytical solution for this problem:
$$ y(t) = y_0 + v_0 t - \frac{g}{2}t^2 $$Let's solve this problem numerically and compare the results.
A second-order ODE can be transformed into two first-order ODE, introducing a new variable:
$$ \dot{y} = v $$$$ \dot{v} = a $$And rewriting Newton's second law as a couple of equations:
$$ \left\{ \begin{array}{r} \frac{\mathrm{d} y}{\mathrm{d}t} = &v, \quad y(t_0) = y_0 \\ \frac{\mathrm{d} v}{\mathrm{d}t} = &-g, \quad v(t_0) = v_0 \end{array} \right.$$First, let's import the necessary Python libraries and customize the environment:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['lines.linewidth'] = 3
matplotlib.rcParams['font.size'] = 13
matplotlib.rcParams['lines.markersize'] = 5
matplotlib.rc('axes', grid=False, labelsize=14, titlesize=16, ymargin=0.05)
matplotlib.rc('legend', numpoints=1, fontsize=11)
This is the equation for calculating the ball trajectory given the model and using the Euler method:
def ball_euler(t0, tend, y0, v0, h):
t, y, v, i = [t0], [y0], [v0], 0
a = -9.8
while t[-1] <= tend and y[-1] > 0:
y.append(y[-1] + h*v[-1])
v.append(v[-1] + h*a)
i += 1
t.append(i*h)
return np.array(t), np.array(y), np.array(v)
Initial values:
y0 = 1
v0 = 20
a = -9.8
Let's call the function with different step sizes:
t100, y100, v100 = ball_euler(0, 10, y0, v0, 0.1)
t10, y10, v10 = ball_euler(0, 10, y0, v0, 0.01)
Here are the plots for the results:
def plots(t100, y100, v100, t10, y10, v10, title):
"""Plots of numerical integration results.
"""
a = -9.8
fig, axs = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(10, 5))
axs[0, 0].plot(t10, y0 + v0*t10 + 0.5*a*t10**2, color=[0, 0, 1, .7], label='Analytical')
axs[0, 0].plot(t100, y100, '--', color=[0, 1, 0, .7], label='h = 100ms')
axs[0, 0].plot(t10, y10, ':', color=[1, 0, 0, .7], label='h = 10ms')
axs[0, 1].plot(t10, v0 + a*t10, color=[0, 0, 1, .5], label='Analytical')
axs[0, 1].plot(t100, v100, '--', color=[0, 1, 0, .7], label='h = 100ms')
axs[0, 1].plot(t10, v10, ':', color=[1, 0, 0, .7], label='h = 10ms')
axs[1, 0].plot(t10, y0 + v0*t10 + 0.5*a*t10**2 - (y0 + v0*t10 + 0.5*a*t10**2),
color=[0, 0, 1, .7], label='Analytical')
axs[1, 0].plot(t100, y100 - (y0 + v0*t100 + 0.5*a*t100**2), '--',
color=[0, 1, 0, .7], label='h = 100ms')
axs[1, 0].plot(t10, y10 - (y0 + v0*t10 + 0.5*a*t10**2), ':',
color=[1, 0, 0, .7], label='h = 10ms')
axs[1, 1].plot(t10, v0 + a*t10 - (v0 + a*t10), color=[0, 0, 1, .7], label='Analytical')
axs[1, 1].plot(t100, v100 - (v0 + a*t100), '--', color=[0, 1, 0, .7], label='h = 100ms')
axs[1, 1].plot(t10, v10 - (v0 + a*t10), ':', color=[1, 0, 0, .7], label='h = 10ms')
ylabel = ['y [m]', 'v [m/s]', 'y error [m]', 'v error [m/s]']
axs[0, 0].set_xlim(t10[0], t10[-1])
axs[1, 0].set_xlabel('Time [s]')
axs[1, 1].set_xlabel('Time [s]')
axs[0, 1].legend()
axs = axs.flatten()
for i, ax in enumerate(axs):
ax.set_ylabel(ylabel[i])
plt.suptitle('Kinematics of a soccer ball - %s method'%title, y=1.02, fontsize=16)
plt.tight_layout()
plt.show()
plots(t100, y100, v100, t10, y10, v10, 'Euler')
Let's use the integrator lsoda
to solve the same problem:
from scipy.integrate import odeint, ode
def ball_eq(yv, t):
y = yv[0] # position
v = yv[1] # velocity
a = -9.8 # acceleration
return [v, a]
yv0 = [1, 20]
t10 = np.arange(0, 4, 0.1)
yv10 = odeint(ball_eq, yv0, t10)
y10, v10 = yv10[:, 0], yv10[:, 1]
t100 = np.arange(0, 4, 0.01)
yv100 = odeint(ball_eq, yv0, t100)
y100, v100 = yv100[:, 0], yv100[:, 1]
plots(t100, y100, v100, t10, y10, v10, 'lsoda')
Let's use an explicit runge-kutta method of order (4)5 due to Dormand and Prince (a.k.a. ode45 in Matlab):
def ball_eq(t, yv):
y = yv[0] # position
v = yv[1] # velocity
a = -9.8 # acceleration
return [v, a]
def ball_sol(fun, t0, tend, yv0, h):
f = ode(fun).set_integrator('dopri5')
# or f = ode(fun).set_integrator('dopri5', nsteps=1, max_step=h/2)
f.set_initial_value(yv0, t0)
data = []
while f.successful() and f.t < tend:
f.integrate(f.t + h)
# or f.integrate(tend)
data.append([f.t, f.y[0], f.y[1]])
data = np.array(data)
return data
data = ball_sol(ball_eq, 0, 4, [1, 20], .1)
t100, y100, v100 = data[:, 0], data[:, 1], data[:, 2]
data = ball_sol(ball_eq, 0, 4, [1, 20], .01)
t10, y10, v10 = data[:, 0], data[:, 1], data[:, 2]
plots(t100, y100, v100, t10, y10, v10, 'dopri5 (ode45)')
Let's consider the air resistance in the calculations for the vertical trajectory of the football ball.
According to the Laws of the Game from FIFA, the ball is spherical, has a circumference of $0.69m$, and a mass of $0.43kg$.
We will model the magnitude of the drag force due to the air resistance by:
Where $\rho$ is the air density $(1.22kg/m^3)$, $A$ the ball cross sectional area $(0.0379m^2)$, and $C_d$ the drag coefficient, which for now we will consider constant and equal to $0.25$ (Bray and Kerwin, 2003).
Applying Newton's second law of motion to the new problem:
In the equation above, $-v/||v||$ takes into account that the drag force always acts opposite to the direction of motion.
Reformulating the second-order ODE above as a couple of first-order equations:
Although (much) more complicated, it's still possible to find an analytical solution for this problem. But for now let's explore the power of numerical integration and use the lsoda
method (the most simple method to call in terms of number of lines of code) to solve this problem:
def ball_eq(yv, t):
g = 9.8 # m/s2
m = 0.43 # kg
rho = 1.22 # kg/m3
cd = 0.25 # dimensionless
A = 0.0379 # m2
y = yv[0] # position
v = yv[1] # velocity
a = -g - 1/(2*m)*rho*cd*A*v*np.abs(v) # acceleration
return [v, a]
yv0 = [1, 20]
t10 = np.arange(0, 4, 0.01)
yv10 = odeint(ball_eq, yv0, t10)
y10, v10 = yv10[:, 0], yv10[:, 1]
def plots(t10, y10, v10):
"""Plots of numerical integration results.
"""
a = -9.8
fig, axs = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(10, 5))
axs[0, 0].plot(t10, y0 + v0*t10 + 0.5*a*t10**2, color=[0, 0, 1, .7], label='No resistance')
axs[0, 0].plot(t10, y10, '-', color=[1, 0, 0, .7], label='With resistance')
axs[0, 1].plot(t10, v0 + a*t10, color=[0, 0, 1, .7], label='No resistance')
axs[0, 1].plot(t10, v10, '-', color=[1, 0, 0, .7], label='With resistance')
axs[1, 0].plot(t10, y0 + v0*t10 + 0.5*a*t10**2 - (y0 + v0*t10 + 0.5*a*t10**2),
color=[0, 0, 1, .7], label='Real')
axs[1, 0].plot(t10, y10 - (y0 + v0*t10 + 0.5*a*t10**2), '-',
color=[1, 0, 0, .7], label='h=10 ms')
axs[1, 1].plot(t10, v0 + a*t10 - (v0 + a*t10), color=[0, 0, 1, .7], label='No resistance')
axs[1, 1].plot(t10, v10 - (v0 + a*t10), '-', color=[1, 0, 0, .7], label='With resistance')
ylabel = ['y [m]', 'v [m/s]', 'y diff [m]', 'v diff [m/s]']
axs[1, 0].set_xlabel('Time [s]')
axs[1, 1].set_xlabel('Time [s]')
axs[0, 1].legend()
axs = axs.flatten()
for i, ax in enumerate(axs):
ax.set_ylabel(ylabel[i])
plt.suptitle('Kinematics of a soccer ball - effect of air resistance', y=1.02, fontsize=16)
plt.tight_layout()
plt.show()
plots(t10, y10, v10)