Here a first-order ODE \(y'=f(t,y)\) is defined over the domain \(t\in [t_0 , t_{\max}]\) and the initial solution \(y(t_0)\) is the known value \(y_0\).
When we speak of solving and ODE we usually mean solving an IVP.
The Euler Method
The Euler method for solving the initial value problem \(y' = f(t, y)\), \(t \in [t_0, t_{\max}]\), \(y(t_0) = y_0\) is
\[y_{n+1} = y_n + h f(t_n ,y_n),\]
where \(h = t_{n+1} - t_n\).
Derivation of the Euler Method
Consider the general first-order IVP
\[ y' = f(t, y), \qquad y(t_0) = y_0, \]
where the solution is \(y(t)\). The tangent line for \(y(t)\) where \(t = t_0\) has the gradient \(f(t_0, y_0)\).
The equation of this tangent line is
\[ y = y_0 + f(t_0, y_0) (t - t_0). \]
Approximating \(y(t_1)\) by calculating the point on the tangent where \(t = t_1\)
\[ y_1 = y_0 + f(t_0, y_0) (t_1 - t_0). \]
We can approximate \(y(t_2)\) using the the tangent line at \(t = t_1\) which has gradient of this \(f(t_1, y_1)\)
\[ y_2 = y_1 + f(t_1, y_1) (t_2 - t_1). \]
Doing similar for more points along the \(t\) domain
So the solution to this IVP using the Euler method is
\(t_n\)
\(y_n\)
0.00
1.000000
0.20
1.000000
0.40
1.040000
0.60
1.123200
0.80
1.257984
1.00
1.459261
Coding the Euler Method (Python)
Import NumPy and matplotlib libraries
import numpy as npimport matplotlib.pyplot as plt
Define Euler method function to solve an IVP
def euler(f, tspan, y0, h): nsteps =int((tspan[1] - tspan[0]) / h) t = np.zeros(nsteps +1) y = np.zeros(nsteps +1) t[0] = tspan[0] y[0] = y0for n inrange(nsteps): y[n+1] = y[n] + h * f(t[n], y[n]) t[n+1] = t[n] + hreturn t, y
Setup IVP and solve using the euler() function.
# Define the ODE functiondef f(t, y):return t * y# Define IVP parameterstspan = [0, 1] # boundaries of the t domainy0 = [1] # initial value of the solutionh =0.2# step length# Solve using the Euler methodt, y = euler(f, tspan, y0, h)# Print table of solution valuesprint("| t | y | ")print("|:----:|:---------:|")for n inrange(len(t)):print(f"| {t[n]:4.2f} | {y[n]:9.6f} |")# Plot solutionfig, ax = plt.subplots()plt.plot(t, y, "bo-", label="Euler")plt.xlabel("$t$")plt.ylabel("$y(t)$")plt.show()
Functions must be defined at the bottom of the script file in MATLAB.
Setup IVP and solve using the euler() function.
% Define ODE functionf=@(t,y) t*y;% Define IVP parameterstspan= [0,1];% boundaries of the t domainy0=1;% initial value of the solutionh=0.2;% step length% Solve IVP using the Euler method[t,y] =euler(f,tspan,y0,h);% Print table of solution values (for loop is used to group print statements)fori=1:1fprintf("| t | y |\n|:----:|:---------:|");forn=1:length(t)fprintf("\n| %4.2f | %9.6f | %9.6f | %8.2e |",t(n),y(n));endend% Plot solutionplot(t,y,"b-o",LineWidth=2,MarkerFaceColor="b")axispaddedxlabel("$t$",FontSize=14,Interpreter="latex")ylabel("$y(t)$",FontSize=14,Interpreter="latex")
where \(\tilde{t_n}\) is some point between \(t_{n}\) and \(t_{n} + h\). Since \(y'(t_n) = f(t_n, y_n)\), substituting equation (Equation 2) into equation (Equation 1) gives
So the local truncation error for the Euler method is proportional to \(h^2\).
The actual value of the local truncation error will change depending on the solution to the ODE and the value of \(t\).
Assuming that the functions \(f\) and its partial derivatives \(f_t\) and \(f_y\) are bounded, we can introduce a constant \(M\) such that
\[ | e_n | \leq \frac{Mh^2}{2}. \]
This means that the worst possible truncation error for the Euler method is \(\dfrac{M h^2}{2}\).
Global Truncation Error (GTE)
The global truncation error denoted by \(E_n\), is the accumulation of the local truncation errors for the steps of the method to compute that solution from \(t=t_0\) up to \(t = t_n\)
We saw that the local truncation error at each step is at most \(\frac{Mh^2}{2}\) then after \(n\) steps the upper bound of the global truncation error is at most \(\frac{nMh^2}{2}\).
Using a constant step length \(h\) then
\[n = \frac{t_n - t_0}{h}, \]
so the upper bound of the global truncation error is
If \(C = \frac{(t_n - t_0 ) M}{2}\) is some positive constant then \(E_n \leq C h\) so we say that \(E_n = O(h)\).
So as the step length \(h\) decreases, the GTE will converge to zero in a linear fashion.
This confirms what we saw in the plot of the errors against step length for the Euler method (Figure 1) where the errors formed a line with a gradient of about 1.
Systems of ODEs
A system of \(N\) first-order ODEs is expressed as a set of multiple ODEs such that
where \(y_1(t), y_2(t), \ldots, y_N(t)\) are multiple dependent functions.
An IVP involving a system of ODEs requires an initial value for each equation in the system, i.e., \(y_1(t_0) = a_1\), \(y_2(t_0) = a_2\) etc. where \(a_1, a_2, \ldots, a_N\) are constants.
Solving a System of ODEs using the Euler Method
To apply a the Euler to solve a system of ODEs, we can write it in vector form \(\mathbf{y}' = f(t, \mathbf{y})\) where
The SIR model is a simple model that describes the spread of an infectious disease in a population.
The model divides the population into three compartments: those who have not yet contracted the disease but are susceptible (\(S\)) those who are infected (\(I\)) and those who have recovered and gained immunity (\(R\)).
The formulation of the model is
\[ \begin{align*}
\frac{\mathrm{d}y}{\mathrm{d}t} &= -\dfrac{\beta}{N} IS, \\
\frac{\mathrm{d}y}{\mathrm{d}t} & =\dfrac{\beta}{N} IS - \gamma I, \\
\frac{\mathrm{d}y}{\mathrm{d}t} &= \gamma I,
\end{align*} \]
where \(N = S + I + R\) is the total population, \(\beta\) is the infection rate in the number of people per day who become infected, and \(\gamma\) is the recovery rate at which a person who is infected recovers per day.
A disease breaks out in a population where \(\beta = 0.5\), \(\gamma = 0.1\) and \(S(0) = 99\), \(I(0) = 1\) and \(R(0) = 0\). Calculate the solution to the SIR model over the first 50 days using the Euler method with a step length of \(h = 1\).
Solution
Let \(y_1 = S\), \(y_2 = I\) and \(y_3 = R\) and writing the SIR model in vector form we have
\[ \begin{align*}
\mathbf{y} &= \begin{pmatrix} S \\ I \\ R \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix}, &
f(t, \mathbf{y}) &=
\begin{pmatrix}
-\dfrac{\beta}{N} IS \\
\dfrac{\beta}{N} IS - \gamma I \\
\gamma I
\end{pmatrix} =
\begin{pmatrix}
-\dfrac{\beta}{N} y_1 y_2 \\
\dfrac{\beta}{N} y_1 y_2 - \gamma y_2 \\
\gamma y_3
\end{pmatrix}.
\end{align*} \]
The initial conditions are \(\mathbf{y}_0 = (99, 1, 0)\) and \(N = 100\), \(\beta = 0.5\), \(\gamma = 0.1\) and \(\frac{\beta}{N} = 0.05\).
The plot of the solution to the SIR model over the first 50 days is shown in Figure 2.
Figure 2: Plot of the solution to the SIR model using the Euler method.
Code (Python)
def euler(f, tspan, y0, h): nsteps =int((tspan[1] - tspan[0]) / h) N =len(y0) t = np.zeros(nsteps +1) y = np.zeros((N, nsteps +1)) t[0] = tspan[0] y[:,0] = y0for n inrange(nsteps): y[:,n+1] = y[:,n] + h * f(t[n], y[:,n]) t[n+1] = t[n] + hreturn t, y.T
Define the SIR model function, setup the IVP and solve using the Euler method.
# Define SIR modeldef SIR(_, y): S, I, R = y N = S + I + R dS =-beta / N * I * S dI = beta / N * I * S - gamma * I dR = gamma * Ireturn np.array([ dS, dI, dR ])# Define IVP parameterstspan = [0, 50] # boundaries of the t domainy0 = [99, 1, 0] # initial valuesbeta, gamma =0.5, 0.1# Model parametersh =1# step length# Solve the IVP using the Euler methodt, y = euler(SIR, tspan, y0, h)# Plot solutionfig, ax = plt.subplots()plt.plot(t, y[:,0], "b-", label="Susceptible")plt.plot(t, y[:,1], "r-", label="Infected")plt.plot(t, y[:,2], "g-", label="Recovered")plt.xlabel("time (days)")plt.ylabel("Population")plt.legend()plt.show()
Define the SIR model function, setup the IVP and solve using the Euler method.
% Define IVP parameterstspan= [0,50];% boundaries of the t domainy0= [99,1,0]';% initial values [S, I, R]beta=0.5;% infection rategamma=0.1;% recovery rateh=1;% step length% Solve the IVP using the Euler method[t,y] =euler(@(t,y)SIR(t,y,beta,gamma),tspan,y0,h);% Plot solutionplot(t,y(:,1),'b-',LineWidth=2)holdonplot(t,y(:,2),'r-',LineWidth=2)plot(t,y(:,3),'g-',LineWidth=2)holdoffxlabel('$t$',Fontsize=14,Interpreter='latex')ylabel('$y$',Fontsize=14,Interpreter='latex')legend('Susceptible','Infected','Recovered',Location="east")axispadded
Higher-order ODEs
The Euler method can only be used to solve first-order ODEs.
For higher order ODEs we rewrite them as system of first-order ODEs and apply the Euler method to solve the system.
For example consider the \(N\)-th order ODE
\[ y^{(N)} = f(t, y, y', y'' ,\ldots ,y^{(N-1)}). \]
If we introduce functions \(y_1(t), y_2(t), \ldots, y_N(t)\) where \(y_1=y\), \(y_2 =y'\), \(y_3 =y''\) and so on then
A mass-spring-damper model consists of objects connected via springs and dampers. A simple example of the application of a model is a single object connected to a surface is shown in Figure 3.
Figure 3: The mass-spring-damper model.
The displacement of the object over time can be modelled by the second-order ODE
\[ m \ddot{y} + c \dot{y} + k y = 0, \]
where \(y\) represents the displacement of the object, \(\dot{y}\) represents the change in the displacement over time (i.e., velocity), \(m\) represents the mass of the object, \(c\) is the damping coefficient and \(k\) is a spring constant based on the length and elasticity of the spring.
An object of mass 1 kg is connected to a dampened spring with \(c = 2\) and \(k = 4\). The object is displaced by 1m and then released.
Use the Euler method to compute the displacement of the object over the first 5 seconds after it was released.
Solution
Let \(y_1 = y\) and \(y_2 = \dot{y}\) then we have the system of two first-order ODEs
\[ \begin{align*}
\dot{y}_1 &= y_2, \\
\dot{y}_2 &= \frac{1}{m}(- c y_2 - k y_1).
\end{align*} \]
Since \(m = 1\), \(c = 2\) and \(k = 4\) then writing the system in vector form \(\dot{\mathbf{y}} = f(t, \mathbf{y})\) gives
