Engineering at Alberta Courses » Backward (Implicit) Euler Method

Solution Methods for IVPs: Backward (Implicit) Euler Method

Backward (Implicit) Euler Method

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 relate the value of $x$ at $t=t_{i+1}$ , namely $x(t_{i+1})$ with $h=t_{i+1}-t_i$ . However, unlike the explicit Euler method, we will use the Taylor series around the point $x(t_{i+1})$ , that is:

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

Substituting the differential equation into the above equation yields:

$x_i=x(t_{i+1})-F(x_{i+1},t_{i+1}) (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+1},t_{i+1}) (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. The backward Euler method is termed an “implicit” method because it uses the slope $\frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t)$ at the unknown point $x_{i+1}$ , namely: $F(x_{i+1},t_{i+1})$ . The developed equation can be linear in $x_{i+1}$ or nonlinear. Nonlinear equations can often be solved using the fixed-point iteration method or the Newton-Raphson method to find the value of $x_{i+1}$ . While the implicit scheme does not provide better accuracy than the explicit scheme, it comes with additional computations. However, one advantage is that this scheme is always stable.

The following Mathematica code adopts the implicit Euler scheme and uses the built-in FindRoot function to solve for $x_{i+1}$ . The code is similar to the code provided for the explicit scheme except for the line to calculate $x_{i+1}$ .

View Mathematica Code

ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
 xtable = Table[0, {i, 1, n}];
 xtable[[1]] = x0;
 Do[  s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}];  xtable[[i]] = x /. s[[1]],  {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
ImplicitEulerMethod[function1, 35, 1.0, 0, 50] // MatrixForm
ImplicitEulerMethod[function2, 35, 1.0, 0, 50] // MatrixForm

View Python Code

import numpy as np
import sympy as sp
x = sp.symbols('x')
def ImplicitEulerMethod(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):
    s = sp.nsolve(x - xtable[i - 1] - h*fp(x, t0 + i*h), x, xtable[i - 1])  
    xtable[i] = s
  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
print(np.vstack(ImplicitEulerMethod(function1, 35, 1.0, 0, 50)))
print(np.vstack(ImplicitEulerMethod(function2, 35, 1.0, 0, 50)))

Example 1

Repeat Example 1 above using the implicit Euler method.

Solution

When $k=0.015$ , the IVP is given by:

$\frac{\mathrm{d}x}{\mathrm{d}t}=0.015 x$

Given a value for $x_i$ at $t_i$ , and taking $h=t_{i+1}-t_i=1$ , the estimate for $x_{i+1}$ can be calculated according to the formula:

$x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(0.015x_{i+1})$

Rearranging yields:

$x_{i+1}=\frac{x_i}{1-0.015}$

With $x_0=35$ million and $t_0=0$ years, the estimate for the population at $t=1$ years using the implicit Euler method is given by:

$x_1=\frac{35}{1-0.015}=35.533$

at $t=2$ years, the value of $x_2$ in millions is given by:

$x_2=\frac{35.533}{1-0.015}=36.074$

Proceeding iteratively, the following is the Excel table showing the values of the population up to $t=50$ years.

When $k=0.015$ , the implicit method provides an error similar to that of the explicit method:

$\begin{split}E&=35e^{(0.015\times 50)}-74.517=74.095-74.517=-0.4\\E_r&=\frac{-0.4}{74.095}=-0.005\end{split}$

When $k=0.05$ , the IVP is given by:

$\frac{\mathrm{d}x}{\mathrm{d}t}=0.05 x$

Given a value for $x_i$ at $t_i$ , and taking $h=t_{i+1}-t_i=1$ , the estimate for $x_{i+1}$ can be calculated according to the formula:

$x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(0.05x_{i+1})$

Rearranging yields:

$x_{i+1}=\frac{x_i}{1-0.05}$

With $x_0=35$ million and $t_0=0$ years, the estimate for the population at $t=1$ years using the implicit Euler method is given by:

$x_1=\frac{35}{1-0.05}=36.842$

at $t=2$ years, the value of $x_2$ in millions is given by:

$x_2=\frac{36.842}{1-0.05}=38.781$

Proceeding iteratively, the following is the Excel table showing the values of the population up to $t=50$ years.

Again, the error produced by the implicit method is comparable to that produced by the explicit method shown above:

$\begin{split}E&=35e^{(0.05\times 50)}-454.871=426.387-454.871=-28.483\\E_r&=\frac{-28.483}{426.387}=-0.07\end{split}$

The following is the graph showing the exact solution versus the data using both the implicit (blue dots) and the explicit (black dots) Euler method for both values of $k$ .

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)

ImplicitEulerMethod[fp_, x0_, h_, t0_,  tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

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];
Data2 = ImplicitEulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data], Blue, Point[Data2]}, 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];
Data2 = ImplicitEulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data], Blue, Point[Data2]}, 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

def ImplicitEulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  x = sp.symbols('x')
  for i in range(1,n):
    s = sp.nsolve(x - xtable[i - 1] - h*fp(x, t0 + i*h), x, xtable[i - 1])  
    xtable[i] = s
  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)
Data2 = ImplicitEulerMethod(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.scatter([point[0] for point in Data2],[point[1] for point in Data2])
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)
Data2 = ImplicitEulerMethod(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.scatter([point[0] for point in Data2],[point[1] for point in Data2])
plt.xlabel("Time (years)"); plt.ylabel("Population (Millions)")
plt.grid(); plt.show()

MATLAB files: File 1 (ex12_4.m) File 2 (implicit_euler.m)

The following tool illustrates the effect of choosing the step size $h$ on the difference between the numerical solution obtained using the implicit methods (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

Repeat Example 2 above using the implicit Euler method.

Solution

Recall that the IVP for this problem is given by:

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

Taking $h=1$ day, $c_{i+1}$ can be predicted for a given value of $c_i$ using the implicit Euler method as follows:

$c_{i+1}=c_i+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_{i+1},t=t_{i+1}}=c_i+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}\left(5e^{\left(\frac{-2t_{i+1}}{1000}\right)}-c_{i+1}\right)$

Rearranging:

$c_{i+1}=\frac{c_i+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}\left(5e^{\left(\frac{-2t_{i+1}}{1000}\right)}\right)}{1+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}}$

With $t_0=0$ and $c_0=5$ ppb, $c_1$ is calculated as:

$c_1=\frac{5+0.015\left(5e^{\left(\frac{-2(1)}{1000}\right)}\right)}{1+0.015}=4.9999$

Proceeding iteratively produces the following graph with results almost identical to those produced by the explicit method. You can download the Microsoft Excel file example2implicit.xlsx showing the calculations.

View Mathematica Code

ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {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 = ImplicitEulerMethod[fp, 5, 1, 0, 365*2];
ListPlot[Data, AxesLabel -> {"time (days)", "Concentration ppb"}]

View Python Code

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

def ImplicitEulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  x = sp.symbols('x')
  for i in range(1,n):
    s = sp.nsolve(x - xtable[i - 1] - h*fp(x, t0 + i*h), x, xtable[i - 1])  
    xtable[i] = s
  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 = ImplicitEulerMethod(fp, 5, 1, 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()

MATLAB files: File 1 (ex12_5.m) File 2 (implicit_euler.m)

Example 3

We will repeat Example 3 above to illustrate that the implicit Euler method is always stable.
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

Using an implicit 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}}=x_i-h x_{i+1}\Rightarrow x_{i+1}=\frac{x_i}{(1+h)}=\frac{x_{i-1}}{(1+h)(1+h)}$

Therefore:

$x_{i+1}=\frac{x_0}{(1+h)^{i+1}}$

$h$ is a positive number. Therefore, the term $\frac{1}{(1+h)^{i+1}}$ is bounded for every $i$ and therefore, $x_{i+1}$ is bounded. The following tool illustrates the difference between the obtained result for $t\in[0,10]$ . Vary the value of $h$ and notice how the obtained solution (black dots) compares to the exact solution (blue curve). While the accuracy is lost with higher values of $h$ , the solution maintains numerical stability.

Example 4

Consider the following nonlinear IVP:

$\frac{\mathrm{d}x}{\mathrm{d}t}=tx^2+2x$

with the initial condition $x(0)=-5$ . Find the values of $x$ for $t\in[0,5]$ taking $h=0.1$ using the implicit Euler method.

Solution

Using Mathematica, the exact solution to the differential equation given can be obtained as:

$x=\frac{-20e^{2t}}{9-5e^{2t}+10e^{2t}t}$

The implicit Euler scheme provides the following estimate for $x_{i+1}$ :

$x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(t_{i+1}x_{i+1}^2+2x_{i+1})$

Since $x_{i+1}$ appears on both sides, and the equation is nonlinear in $x_{i+1}$ , therefore, the Newton-Raphson method will be used at each time increment to find the value of $x_{i+1}$ !
Setting $x_0=-5$ , $h=0.1$ , $t_0=0$ , and $t_1=0.1$ yields:

$x_1=-5+0.1(0.1 x_1^2+2x_1)$

Solving the above nonlinear equation in $x_1$ using the Newton-Raphson method yields:

$x_1=-5.82576$

Similarly, with $t_2=0.2$ , the equation for $x_2$ is:

$x_2=-5.82576+0.1(0.2x_2^2+2x_2)$

Using the Newton-Raphson method:

$x_2=-6.29235$

Proceeding iteratively produces the following table:

The following is the plot of the exact solution (blue line) with the data obtained using the implicit Euler method. The Mathematica code is given below.

View Mathematica Code

ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = ImplicitEulerMethod[fp, -5.0, 0.1, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]},AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm

View Python Code

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
sp.init_printing(use_latex=True)

def ImplicitEulerMethod(fp, x0, h, t0, tmax):
  n = int((tmax - t0)/h + 1)
  xtable = [0 for i in range(n)]
  xtable[0] = x0
  x = sp.symbols('x')
  for i in range(1,n):
    s = sp.nsolve(x - xtable[i - 1] - h*fp(x, t0 + i*h), x, xtable[i - 1])  
    xtable[i] = s
  Data = [[t0 + i*h, xtable[i]] for i in range(n)]
  return Data

x = sp.Function('x')
t = sp.symbols('t')
sol = sp.dsolve(x(t).diff(t) - t*x(t)**2 - 2*x(t), ics={x(0): -5})
display(sol)
def fp(x, t): return t*x**2 + 2*x
Data = ImplicitEulerMethod(fp, -5.0, 0.1, 0, 5)
x_val = np.arange(0,5,0.01)
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],c='k')
plt.xlabel("t"); plt.ylabel("x")
plt.grid(); plt.show()
print(["t_i", "x_i"],"\n",np.vstack(Data))

MATLAB files: File 1 (ex12_7.m) File 2 (implicit_euler.m)