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
![Rendered by QuickLaTeX.com x](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-822065bc13c102457e5826bb62632b02_l3.png)
![Rendered by QuickLaTeX.com x_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-6e607e184ef7ff578cf200ef43f3f987_l3.png)
![Rendered by QuickLaTeX.com t_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-e21115ab3d9ad47abb0226449b97426f_l3.png)
![Rendered by QuickLaTeX.com x](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-822065bc13c102457e5826bb62632b02_l3.png)
![Rendered by QuickLaTeX.com t=t_{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-34bd47dd228c2048f61284ee53d1e0ef_l3.png)
![Rendered by QuickLaTeX.com x(t_{i+1})](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-bb066f27a30f214b986a6c18487fdee4_l3.png)
![Rendered by QuickLaTeX.com h=t_{i+1}-t_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-2f576dca3394f4e7e11482c159349565_l3.png)
Substituting the differential equation into the above equation yields:
Therefore, as an approximation, an estimate for
![Rendered by QuickLaTeX.com x(t_{i+1})](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-bb066f27a30f214b986a6c18487fdee4_l3.png)
![Rendered by QuickLaTeX.com x_{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-242ebad07a9a05752788317b0eb627aa_l3.png)
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
![Rendered by QuickLaTeX.com x](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-822065bc13c102457e5826bb62632b02_l3.png)
If the errors from each interval are added together, with
![Rendered by QuickLaTeX.com n=\frac{L}{h}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-de81fd39e2c0320c1fa6fa6aef07f745_l3.png)
![Rendered by QuickLaTeX.com L](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-8fc4a64bb149669d8ae9594aeec2db6c_l3.png)
![Rendered by QuickLaTeX.com t_n-t_0](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-1f20a20364adda1d0291ac048586af13_l3.png)
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 CodeEulerMethod[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
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
![Rendered by QuickLaTeX.com x(0)=35](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-fd2d96edaa87f81c34310a4dbba04c4c_l3.png)
If the birth rate is 6%, then
![Rendered by QuickLaTeX.com k=0.06-0.01=0.05](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-fd40cbfc6b28ee07b8bf884b247b4b7a_l3.png)
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
![Rendered by QuickLaTeX.com x_2](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-618163b950cd7633ea4bf03c9e3ce090_l3.png)
Iterating further, the values of
![Rendered by QuickLaTeX.com x](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-822065bc13c102457e5826bb62632b02_l3.png)
![Rendered by QuickLaTeX.com t_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-e21115ab3d9ad47abb0226449b97426f_l3.png)
![Rendered by QuickLaTeX.com x_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-6e607e184ef7ff578cf200ef43f3f987_l3.png)
![IVPe1](https://engcourses-uofa.ca/wp-content/uploads/IVPe1.png)
When
![Rendered by QuickLaTeX.com k=0.015](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a995def1d50ea37125a727298310ccac_l3.png)
![Rendered by QuickLaTeX.com t=50](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-f344081b097f2d8e6e01b95576318e77_l3.png)
![Rendered by QuickLaTeX.com t=50](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-f344081b097f2d8e6e01b95576318e77_l3.png)
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
![Rendered by QuickLaTeX.com t_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-e21115ab3d9ad47abb0226449b97426f_l3.png)
![Rendered by QuickLaTeX.com x_i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-6e607e184ef7ff578cf200ef43f3f987_l3.png)
![IVPe2](https://engcourses-uofa.ca/wp-content/uploads/IVPe2.png)
Compared to the previous case, when
![Rendered by QuickLaTeX.com k=0.05](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a39c8714b93c97a01a80a9f60c2a01f8_l3.png)
![Rendered by QuickLaTeX.com t=50](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-f344081b097f2d8e6e01b95576318e77_l3.png)
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
![Rendered by QuickLaTeX.com t](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-4d1ce3eee1a1120ff6d4bb054b3014f0_l3.png)
![IVP2](https://engcourses-uofa.ca/wp-content/uploads/IVP2.png)
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)"}]
# 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 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
![Rendered by QuickLaTeX.com c](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-7855ac73ada79d2b379708df0d5b6fb3_l3.png)
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
![Rendered by QuickLaTeX.com h=t_{i+1}-t_i=1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-f5dcaa117cae08290771df25e7899902_l3.png)
![Rendered by QuickLaTeX.com t_0=0](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-7e6a951d564386d428a51e1874d5374b_l3.png)
![Rendered by QuickLaTeX.com c_0=5](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-79a87f1f0e0131bbddd6a5fb44173238_l3.png)
![Rendered by QuickLaTeX.com c](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-7855ac73ada79d2b379708df0d5b6fb3_l3.png)
![Rendered by QuickLaTeX.com t=730](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-afbeec2224c62611400b12ba5ddb0348_l3.png)
Taking
![Rendered by QuickLaTeX.com h=1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-05e275b0ab6e717d6412c18c3e2ee0a1_l3.png)
![Rendered by QuickLaTeX.com c_1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-19b1488688ebc02786f5a3ad99f15780_l3.png)
![Rendered by QuickLaTeX.com c_2](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-2cf6fbcd2fde4b9e55b8a5dfc025db7c_l3.png)
![Rendered by QuickLaTeX.com t_1=1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-1fafbe928bbb61fc87ac66719c52489e_l3.png)
At
![Rendered by QuickLaTeX.com t_2=2](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-0d3272a7b6e4667cce03253f1dd52bfc_l3.png)
Carrying on produces the values of
![Rendered by QuickLaTeX.com c](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-7855ac73ada79d2b379708df0d5b6fb3_l3.png)
![Rendered by QuickLaTeX.com t=2(365)](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-781cffbe31323ad5461613de13d5ddab_l3.png)
![IPV1](https://engcourses-uofa.ca/wp-content/uploads/IPV1.png)
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"}]
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
![Rendered by QuickLaTeX.com x_0=1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a1821f84ea65f7c9738882d77ed0986e_l3.png)
![Rendered by QuickLaTeX.com h=0.1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-1ca78a11d8ddb68c4b2f3607b6c288c1_l3.png)
![Rendered by QuickLaTeX.com h=3](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a905876bb406a4385242ccf76f7fc3e0_l3.png)
Solution
The exact solution to the given equation is given by:
Using an explicit Euler scheme, the value of
![Rendered by QuickLaTeX.com x_{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-242ebad07a9a05752788317b0eb627aa_l3.png)
Therefore:
When
![Rendered by QuickLaTeX.com h=0.1](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-1ca78a11d8ddb68c4b2f3607b6c288c1_l3.png)
![Rendered by QuickLaTeX.com (1-h)^{i+1}=0.9^{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a11d6d7c9ca72493b39ea8838200f1e7_l3.png)
![Rendered by QuickLaTeX.com i](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-2b061624d8dea4885990bf08e6bd9eed_l3.png)
![Rendered by QuickLaTeX.com x_{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-242ebad07a9a05752788317b0eb627aa_l3.png)
![Rendered by QuickLaTeX.com h=3](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-a905876bb406a4385242ccf76f7fc3e0_l3.png)
![Rendered by QuickLaTeX.com (1-h)^{i+1}=(-2)^{i+1}](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-0b9a73b5457de43d17d6bf710cef32ff_l3.png)
The following tool illustrates the difference between the obtained result for
![Rendered by QuickLaTeX.com t\in[0,10]](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-ff10d62e92b48bda015872665c24c655_l3.png)
![Rendered by QuickLaTeX.com h](https://engcourses-uofa.ca/wp-content/ql-cache/quicklatex.com-bfbf0d614550cb5749f6c2febeda050d_l3.png)
Lecture Video