Open Educational Resources

Solution Methods for IVPs: (Explicit) Euler Method

(Explicit) Euler Method

The Euler method is one of the simplest methods for solving first-order IVPs. Consider the following IVP:

    \[\frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t)\]


Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, we can use a Taylor approximation to estimate the value of x at t=t_{i+1}, namely x(t_{i+1}) with h=t_{i+1}-t_i:

    \[x(t_{i+1})=x_i+\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{t_i} (h) +\mathcal{O}(h^2)\]


Substituting the differential equation into the above equation yields:

    \[x(t_{i+1})=x_i+F(x_i,t_i) (h) +\mathcal{O}(h^2)\]


Therefore, as an approximation, an estimate for x(t_{i+1}) can be taken as x_{i+1} as follows:

    \[x(t_{i+1})\approx x_{i+1}=x_i+F(x_i,t_i) (h)\]


Using this estimate, the local truncation error is thus proportional to the square of the step size with the constant of proportionality related to the second derivative of x, which is the first derivative of the given IVP:

    \[E_{\mbox{local}}=\mathcal{O}(h^2)\]


If the errors from each interval are added together, with n=\frac{L}{h} being the number of intervals and L the total length t_n-t_0, then, the total error is:

    \[E=\mathcal{O}(h^2) \frac{L}{h}=\mathcal{O}(h)\]

There are many examples in engineering and biology in which such IVPs appear. Check this page for examples.

The following Mathematica code provides a procedure whose inputs are

  • The differential equation as a function of the dependent variable x and the independent variable t.
  • The step size h.
  • The initial value t_0 and the final value t_{max} of the independent variable.
  • The initial value x(t_0)=x_0.

The procedure then carries on the Euler method and outputs the required data vector.

View Mathematica Code
EulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] = xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2,  n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)
function1[y_, t_] := 0.05 y
function2[x_, t_] := 0.015 x
EulerMethod[function1, 35, 1.0, 0, 50] // MatrixForm
EulerMethod[function2, 35, 1.0, 0, 50] // MatrixForm
View Python Code
def EulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  for i in range(1,n):
    xtable[i] = xtable[i - 1] + h*fp(xtable[i - 1], t0 + (i - 1)*h)
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data
def function1(y, t): return 0.05*y
def function2(x, t): return 0.015*x
display(EulerMethod(function1, 35, 1.0, 0, 50))
display(EulerMethod(function2, 35, 1.0, 0, 50))

Example 1

The Canadian population at t=0 (current year) is 35 million. If 2.5% of the population have a child in a given year, while the death rate is 1% of the population, what will the population be in 50 years? What if 6% of the population have a child in a given year and the death rate is kept constant at 1%, what will the population be in 50 years?

Solution

The rate of growth is directly proportional to the current population with the constant of proportionality k=0.025-0.01=0.015. If x represents the population in millions and t represents the time in years, then the IVP is given by:

    \[\frac{\mathrm{d}x}{\mathrm{d}t}=0.015 x\]


With the initial condition of x(0)=35, the exact solution to this differential equation is:

    \[x=35e^{0.015t}\]


If the birth rate is 6%, then k=0.06-0.01=0.05 and the exact solution of the population follows the following equation:

    \[x=35e^{0.05t}\]

The Euler method provides a numerical solution as follows. Using a step size of h=t_{i+1}-t_i=1, we can set t_0=0, t_1=1, t_2=2, \cdots, t_{50}=50 years. The value of x_0 is given as 35 million. When k=0.015, the value of x_1 in millions is given by:

    \[x_1=x_0+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_0,t=t_0}=35+0.015\times 35=35.525\]


The value of x_2 is given by:

    \[x_2=x_1+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_1,t=t_1}=35.525+0.015\times 35.525=36.0579\]


Iterating further, the values of x at each time point can be obtained. The following is the Microsoft Excel table showing the values of t_i in years and the corresponding values of the population x_i in millions:
IVPe1
When k=0.015, the numerical procedure produces a very good approximation for the population at t=50 years. The error between the prediction of the exact solution to the differential equation and the numerical value at t=50 is given in millions by:

    \[E=35e^{(0.015\times 50)} - 73.6835=74.095-73.6835=0.4\]


The relative error in this case is given by:

    \[E_r=\frac{0.4}{74.095}=0.005\]

For the second case, when k=0.05, the value of x_1 in millions is given by:

    \[x_1=x_0+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_0,t=t_0}=35+0.05\times 35=36.75\]


Proceeding iteratively, the following is the Microsoft Excel table showing the values of t_i in years and the corresponding values of the population x_i in millions:
IVPe2
Compared to the previous case, when k=0.05, the numerical procedure is not as good. The error between the prediction of the exact solution to the differential equation and the numerical value at t=50 is given in millions by:

    \[E=35e^{(0.05\times 50)} - 401.359=426.387-401.359=25.0283\]


The relative error in this case is given by:

    \[E_r=\frac{25.0283}{426.387}=0.06\]


The following is the graph showing the exact solution versus the data using the Euler method for both cases. The Mathematica code used is given below. Notice how the error in the Euler method increases as t increases.
IVP2

View Mathematica Code
EulerMethod[fp_, x0_, h_, t0_, tmax_] := 
  (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] =  xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

Clear[x, xtable]
k = 0.015;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]

Clear[x, xtable]
k = 0.05;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]
View Python Code
# UPGRADE: need Sympy 1.2 or later, upgrade by running: "!pip install sympy --upgrade" in a code cell
# !pip install sympy --upgrade
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
sp.init_printing(use_latex=True)

def EulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  for i in range(1,n):
    xtable[i] = xtable[i - 1] + h*fp(xtable[i - 1], t0 + (i - 1)*h)
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data

x = sp.Function('x')
t = sp.symbols('t')
k = 0.015
sol = sp.dsolve(x(t).diff(t) - k*x(t), ics={x(0): 35})
display(sol)
def fp(x, t): return k*x
Data = EulerMethod(fp, 35, 1.0, 0, 50)
x_val = np.arange(0,50,0.1)
plt.plot(x_val, [sol.subs(t, i).rhs for i in x_val])
plt.scatter([point[0] for point in Data],[point[1] for point in Data])
plt.xlabel("Time (years)"); plt.ylabel("Population (Millions)")
plt.grid(); plt.show()

k = 0.05
sol = sp.dsolve(x(t).diff(t) - k*x(t), ics={x(0): 35})
display(sol)
Data = EulerMethod(fp, 35, 1.0, 0, 50)
plt.plot(x_val, [sol.subs(t, i).rhs for i in x_val])
plt.scatter([point[0] for point in Data],[point[1] for point in Data])
plt.xlabel("Time (years)"); plt.ylabel("Population (Millions)")
plt.grid(); plt.show()

The following tool illustrates the effect of choosing the step size h on the difference between the numerical solution (shown as dots) and the exact solution (shown as the blue curve) for the case when k=0.05. The smaller the step size, the more accurate the solution is. The higher the value of h, the more the numerical solution deviates from the exact solution. The error increases away from the initial conditions.

     

 

 

 

 

 



Example 2

This example is adopted from this page. A lake of volume 10,000 m^3 has an initial concentration of a particular pollutant of 5 parts per billion (ppb). This pollutant is due to a pesticide that is no longer available in the market. The volume of the lake is constant throughout the year, with daily water flowing into and out of the lake at the rate of q=100+50\cos{(2\pi t/365)} \frac{m^3}{day}. As the pesticide is no longer used, the concentration of the pollutant in the surrounding soil is given in ppb as a function of t(days) by c_{in}=5e^{\left(\frac{-2t}{1000}\right)} which indicates that the concentration is following an exponential decay. This concentration can be assumed to be that of the pollutant in the fluid flowing into the lake. The concentration of the pollutant in the fluid flowing out of the lake is the same as the concentration in the lake. What is the concentration of this particular pollutant in the lake after two years?

Solution

The first step in this problem is to properly write the differential equation that needs to be solved. If \alpha is the total amount of pollutant in the lake, then, the rate of change of \alpha with respect to t is equal to the total amount of pollutant entering into the lake per day. Each day, the total amount of pollutant entering the lake is given by c_{in}\times q, while the total amount of pollutant exiting the lake is given by c\times q, where c is the concentration of the pollutant in the lake. Therefore:

    \[\frac{\mathrm{d}\alpha}{\mathrm{d}t}=\left(100+50\cos{(2\pi t/365)}\right)\left(c_{in}-c\right)\]


The concentration is equal to the total amount divided by the volume of the lake, therefore, the differential equation in terms of c is given by:

    \[\frac{\mathrm{d}c}{\mathrm{d}t}=\frac{\left(100+50\cos{(2\pi t/365)}\right)}{10000}\left(c_{in}-c\right)=\frac{\left(100+50\cos{(2\pi t/365)}\right)}{10000}\left(5e^{\left(\frac{-2t}{1000}\right)}-c\right)\]


Notice that finding an exact solution to the above differential equation is not an easy task and a numerical solution would be the preferred solution in this case. Taking h=t_{i+1}-t_i=1 day, t_0=0, and c_0=5 ppb the Euler method can be used to find the concentration c at t=730 days. You can download the Microsoft Excel file example2.xlsx showing the calculations and the produced graph.
Taking h=1 day, we will show the calculations for the estimates c_1 and c_2. At t_1=1 day:

    \[\begin{split}c_1&=c_0+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_0,t=t_0}\\&=5+1\left(\frac{\left(100+50\cos{(2\pi 0/365)}\right)}{10000}\left(5e^{\left(\frac{-2(0)}{1000}\right)}-5\right)\right)\\&=5\end{split}\]


At t_2=2 days:

    \[\begin{split}c_2&=c_1+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_1,t=t_1}\\&=5+1\left(\frac{\left(100+50\cos{(2\pi 1/365)}\right)}{10000}\left(5e^{\left(\frac{-2(1)}{1000}\right)}-5\right)\right)\\&=4.9999\end{split}\]


Carrying on produces the values of c at the remaining time intervals up to t=2(365) days. The following is the output produced by the Mathematica code below.
IPV1

View Mathematica Code
EulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] =  xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)
f[t_] := 100 + 50*Cos[2 Pi*t/365];
cin[t_] := 5 E^(-2 t/1000);
fp[c_, t_] := f[t] (cin[t] - c)/10000;
Data = EulerMethod[fp, 5, 1.0, 0, 365*2];
ListPlot[Data, AxesLabel -> {"time (days)", "Concentration ppb"}]
View Python Code
import numpy as np
import matplotlib.pyplot as plt

def EulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  for i in range(1,n):
    xtable[i] = xtable[i - 1] + h*fp(xtable[i - 1], t0 + (i - 1)*h)
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data
def f(t): return 100 + 50*np.cos(2*np.pi*t/365)
def cin(t): return 5*np.exp(-2*t/1000)
def fp(c, t): return f(t)*(cin(t) - c)/10000
Data = EulerMethod(fp, 5, 1.0, 0, 365*2)
plt.scatter([point[0] for point in Data],[point[1] for point in Data])
plt.xlabel("Time (days)"); plt.ylabel("Concentration ppb")
plt.grid(); plt.show()

Example 3

In this example we investigate the issue of stability of numerical solution of differential equations using the Euler method. Consider the IVP of the form:

    \[\frac{\mathrm{d}x}{\mathrm{d}t}=-x\]


If x_0=1, compare the two solutions taking a step size of h=0.1, and a step size of h=3.

Solution

The exact solution to the given equation is given by:

    \[x=e^{-t}\]


Using an explicit Euler scheme, the value of x_{i+1} can be obtained as follows:

    \[x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x_i}=(1-h)x_i=(1-h)(1-h)x_{i-1}\]


Therefore:

    \[x_{i+1}=(1-h)^{i+1}x_0\]


When h=0.1, the term (1-h)^{i+1}=0.9^{i+1} is bounded for any value of i which means that the value of x_{i+1} is always bounded. However, if h=3, then, the term (1-h)^{i+1}=(-2)^{i+1} oscillates wildly leading to an instability in the numerical scheme!
The following tool illustrates the difference between the obtained result for t\in[0,10]. Vary the value of h and try to identify the point above which the obtained solution (black dots) starts oscillating wildly around the exact solution (blue curve).

  

 

 

 

 



Lecture Video

 

 


 

 


 

Leave a Reply

Your email address will not be published.