## 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:

Assuming that the value of the dependent variable (say ) is known at an initial value , then, we can use a Taylor approximation to estimate the value of at , namely with :

Substituting the differential equation into the above equation yields:

Therefore, as an approximation, an estimate for can be taken as as follows:

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 , which is the first derivative of the given IVP:

If the errors from each interval are added together, with being the number of intervals and the total length , then, the total error is:

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 and the independent variable .
• The step size .
• The initial value and the final value of the independent variable.
• The initial value .

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 (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 . If represents the population in millions and represents the time in years, then the IVP is given by:

With the initial condition of , the exact solution to this differential equation is:

If the birth rate is 6%, then and the exact solution of the population follows the following equation:

The Euler method provides a numerical solution as follows. Using a step size of , we can set years. The value of is given as 35 million. When , the value of in millions is given by:

The value of is given by:

Iterating further, the values of at each time point can be obtained. The following is the Microsoft Excel table showing the values of in years and the corresponding values of the population in millions:

When , the numerical procedure produces a very good approximation for the population at years. The error between the prediction of the exact solution to the differential equation and the numerical value at is given in millions by:

The relative error in this case is given by:

For the second case, when , the value of in millions is given by:

Proceeding iteratively, the following is the Microsoft Excel table showing the values of in years and the corresponding values of the population in millions:

Compared to the previous case, when , 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 is given in millions by:

The relative error in this case is given by:

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 increases.

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
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 on the difference between the numerical solution (shown as dots) and the exact solution (shown as the blue curve) for the case when . The smaller the step size, the more accurate the solution is. The higher the value of , 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 has an initial concentration of a particular pollutant of 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 . As the pesticide is no longer used, the concentration of the pollutant in the surrounding soil is given in ppb as a function of (days) by 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 is the total amount of pollutant in the lake, then, the rate of change of with respect to 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 , while the total amount of pollutant exiting the lake is given by , where is the concentration of the pollutant in the lake. Therefore:

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

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 day, , and ppb the Euler method can be used to find the concentration at days. You can download the Microsoft Excel file example2.xlsx showing the calculations and the produced graph.
Taking day, we will show the calculations for the estimates and . At day:

At days:

Carrying on produces the values of at the remaining time intervals up to days. The following is the output produced by the Mathematica code below.

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:

If , compare the two solutions taking a step size of , and a step size of .

##### Solution

The exact solution to the given equation is given by:

Using an explicit Euler scheme, the value of can be obtained as follows:

Therefore:

When , the term is bounded for any value of which means that the value of is always bounded. However, if , then, the term oscillates wildly leading to an instability in the numerical scheme!
The following tool illustrates the difference between the obtained result for . Vary the value of and try to identify the point above which the obtained solution (black dots) starts oscillating wildly around the exact solution (blue curve).