Euler’s method - Fundamentals of Numerical Computation

Let a first-order initial-value problem be given in the form

\begin{split} u'(t) &= f\bigl(t,u(t)\bigr), \qquad a \le t \le b,\\ u(a)& =u_0. \end{split}

(6.2.1)

We represent a numerical solution of an IVP by its values at a finite collection of nodes, which for now we require to be equally spaced:

t_i = a + ih, \qquad h=\frac{b-a}{n}, \qquad i=0,\ldots,n.

(6.2.2)

The number $h$ is called the step size.

Because we don’t get exactly correct values of the solution at the nodes, we need to take some care with the notation. From now on we let $\hat{u}(t)$ denote the exact solution of the IVP. The approximate value at $t_i$ computed at the nodes by our numerical methods will be denoted by $u_i\approx \hat{u}(t_i)$ . Because we are given the initial value $u(a)=u_0$ exactly, there is no need to distinguish whether we mean $u_0$ as the exact or the numerical solution.

Consider a piecewise linear interpolant to the (as yet unknown) values $u_0,u_1,\ldots$ , $u_n$ . For $t_i < t < t_{i+1}$ , its slope is

\frac{u_{i+1} - u_{i}}{t_{i+1}-t_i} = \frac{u_{i+1}-u_i}{h}.

(6.2.3)

We can connect this derivative to the differential equation by following the model of $u'=f(t,u)$ :

\frac{u_{i+1}-u_i}{h} = f(t_i,u_i), \qquad i=0,\ldots,n-1.

(6.2.4)

The left-hand side is a forward-difference approximation to $u'(t)$ at $t=t_i$ . We can rearrange the equation to get Euler’s method, our first method for IVPs.

Euler’s method marches ahead in $t$ , obtaining the solution at a new time level explicitly in terms of the latest value.

A basic implementation of Euler’s method is shown in Function 6.2.1.

Algorithm 6.2.1 (euler)

Julia

MATLAB

Python

Euler’s method for an initial-value problem

euler.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
"""
    euler(ivp, n)

Apply Euler's method to solve the given IVP using `n` time steps.
Returns a vector of times and a vector of solution values.
"""
function euler(ivp, n)
    # Time discretization.
    a, b = ivp.tspan
    h = (b - a) / n
    t = [a + i * h for i in 0:n]

    # Initial condition and output setup.
    u = fill(float(ivp.u0), n+1)

    # The time stepping iteration.
    for i in 1:n
        u[i+1] = u[i] + h * ivp.f(u[i], ivp.p, t[i])
    end
    return t, u
end

Euler’s method for an initial-value problem

eulersys.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function [t,u] = eulersys(du_dt, tspan, u0, n)
    % EULERSYS   Euler's method for a first-order IVP system.
    % Input:
    %   du_dt   defines f in u'(t) = f(t, u) 
    %   tspan   endpoints of time interval (2-vector)
    %   u0      initial value (m-vector)
    %   n       number of time steps (integer)
    % Output:
    %   t       selected nodes (vector, length n+1)
    %   u       solution values (array, n+1 by m)

    % Define the time discretization.
    a = tspan(1);  b = tspan(2);
    h = (b - a) / n;
    t = a + (0:n)' * h;

    % Initialize solution array.
    u = zeros(length(u0), n+1);
    u(:, 1) = u0(:);

    % Time stepping.
    for i = 1:n
        u(:, i+1) = u(:, i) + h * du_dt(t(i), u(:, i));
    end

    u = u.';    % conform to MATLAB output convention
end

Euler’s method for an initial-value problem

euler.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def euler(du_dt, tspan, u0, n):
    """
    euler(du_dt, tspan, u0, n)

    Apply Euler's method to solve the IVP u'=du_dt(u,t) over the interval tspan with
    u(tspan[1])=u0, using n subintervals/steps. Return vectors of times and solution
    values.
    """
    a, b = tspan
    h = (b - a) / n
    t = np.linspace(a, b, n+1)
    u = np.zeros((n+1, len(np.atleast_1d(u0))))
    u[0] = u0
    for i in range(n):
        u[i+1] = u[i] + h * du_dt(t[i], u[i])

    return t, u.T

6.2.1Local truncation error¶

Let $\hat{u}(t)$ be the exact solution of the IVP (6.2.1), and suppose that somehow we have access to it at $t=t_i$ , so that $u_i=\hat{u}(t_i)$ . How good is $u_{i+1}$ as an approximation to $\hat{u}(t_{i+1})$ ? The answer is revealed through a Taylor series:

\begin{split} \hat{u}(t_{i+1}) - \bigl[ u_i + hf(t_i,u_i) \bigr] &= \hat{u}(t_{i+1}) - \bigl[ \hat{u}(t_i) + hf\bigl(t_i,\hat{u}(t_i)\bigr) \bigr] \\ &= \bigl[ \hat{u}(t_i) + h \hat{u}'(t_i) + \tfrac{1}{2}h^2 \hat{u}''(t_i) + O(h^3) \bigr] - \bigl[ \hat{u}(t_i) + h\hat{u}'(t_i) \bigr] \notag \\ &= \tfrac{1}{2}h^2 \hat{u}''(t_i) + O(h^3), \end{split}

(6.2.6)

where we used the fact that $\hat{u}$ satisfies the differential equation.

We now introduce some formalities.

Euler’s method is the particular case of (6.2.7) with $\phi(t,u,h) = f(t,u)$ , but we will see other one-step methods in later sections.

In close analogy with Convergence of finite differences, we define truncation error as the residual of (6.2.7) when the exact solution is inserted.

One result follows immediately from the definitions:

6.2.2Convergence¶

While the local truncation error is straightforward to calculate from its definition, it is not the quantity we want to know about and control.

At times the term global error may be interpreted as the max-norm of the global error vector, or as its final value.

By our definitions, the local error in stepping from $t_i$ to $t_{i+1}$ is $h\tau_{i+1}(h)$ . To reach the time $t=b$ from $t=a$ with step size $h$ , we need to take $n=(b-a)/h$ steps. If we want to reach, say, $t=(a+b)/2$ , then we would have to take $n/2$ steps, and so on. In fact, to reach any fixed time in the interval, we need to take $O(n)=O(h^{-1})$ steps. By expressing the local error with a factor of $h$ taken out, the LTE $\tau$ itself is accounting for the simple accumulation of error caused by taking $O(n)$ steps.^[1]

However, global error is not as simple as a sum of local errors. As explained in Theorem 6.1.2 and illustrated in Example 6.1.5, each step causes a perturbation of the solution that can grow as $t$ advances. Thus, we have to account for the flow evolution of individual step truncation errors as well as their mere accumulation. That is the subject of the following theorem.

Proof

Define the global error sequence $ϵ_i=\hat{u}(t_i)-u_i$ . Using (6.2.7), we obtain

ϵ_{i+1} - ϵ_i = \hat{u}(t_{i+1}) - \hat{u}(t_i) - ( {u}_{i+1} - u_i ) = \hat{u}(t_{i+1}) - \hat{u}(t_i) - h\phi(t_i,u_i,h),

(6.2.12)

ϵ_{i+1} = ϵ_i + [\hat{u}(t_{i+1}) - \hat{u}(t_i) - h\phi(t_i,\hat{u}(t_i),h)] + h[\phi(t_i,\hat{u}(t_i),h)- \phi(t_i,u_i,h)].

(6.2.13)

We apply the triangle inequality, (6.2.8), and (6.2.9) to find

|ϵ_{i+1}| \le |ϵ_i| + Ch^{p+1} + h \left| \phi(t_i,\hat{u}(t_i),h)- \phi(t_i,u_i,h)\right|.

(6.2.14)

The Fundamental Theorem of Calculus implies that

\begin{split} \left| \phi(t_i,\hat{u}(t_i),h)- \phi(t_i,u_i,h)\right| & = \left| \int_{u_i}^{\hat{u}(t_i)} \frac{\partial \phi}{\partial u} \,du \right|\\ & \le \int_{u_i}^{\hat{u}(t_i)} \left|\frac{\partial \phi}{\partial u}\right| \,du \\[1mm] & \le L | \hat{u}(t_i)-u_i| = L\, |ϵ_i|. \end{split}

(6.2.15)

Thus

\begin{split} |ϵ_{i+1}| &\le Ch^{p+1} + (1 + hL) |ϵ_i| \\ &\le Ch^{p+1} + (1 + hL) \bigl[ Ch^{p+1} + (1 + hL) |ϵ_{i-1}| \bigr]\\ &\;\vdots \\ &\le Ch^{p+1} \left[ 1 + (1+hL) + (1+hL)^2 + \cdots + (1+hL)^i \right]. \end{split}

(6.2.16)

To get the last line we applied the inequality recursively until reaching $ϵ_0$ , which is zero. Replacing $i+1$ by $i$ and simplifying the geometric sum, we get

|ϵ_i| \le Ch^{p+1}\frac{(1+hL)^i - 1}{(1+hL)-1} = \frac{Ch^p}{L} \left[ (1+hL)^i - 1 \right].

(6.2.17)

We observe that $1+x \le e^x$ for $x\ge 0$ (see Exercise 6.2.5). Hence $(1+hL)^i \le e^{i h L}$ , which completes the proof.

The theorem justifies one more general definition.

We could restate Theorem 6.2.1 as saying that the global error has the same order of accuracy as the LTE. Note, however, that the $O(h^p)$ convergence hides a leading constant that grows exponentially in time. When the time interval is bounded as $h\to 0$ , this does not interfere with the conclusion, but the behavior as $t\to\infty$ contains no such guarantee.

Example 6.2.1 (Convergence of Euler’s method)

Julia

MATLAB

Python

Example 6.2.1

We consider the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

using OrdinaryDiffEq
f(u, p, t) = sin((t + u)^2);
tspan = (0.0, 4.0);
u0 = -1.0;
ivp = ODEProblem(f, u0, tspan)

ODEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 4.0)
u0: -1.0

Here is the call to Function 6.2.1.

using Plots
t, u = FNC.euler(ivp, 20)
plot(t, u;
    m=2,  label="n=20", 
    xlabel=L"t",  ylabel=L"u(t)",
    title="Solution by Euler's method")

We could define a different interpolant to get a smoother picture above, but the derivation of Euler’s method assumed a piecewise linear interpolant, and that sets the limit of its accuracy. We can instead request more steps to make the interpolant look smoother.

t, u = FNC.euler(ivp, 50)
plot!(t, u, m=2, label="n=50")

Increasing $n$ changed the solution noticeably. Since we know that interpolants and finite differences become more accurate as $h\to 0$ , we should anticipate the same behavior from Euler’s method. We don’t have an exact solution to compare to, so we will use a DifferentialEquations solver to construct an accurate reference solution.

u_exact = solve(ivp, Tsit5(), reltol=1e-14, abstol=1e-14)
plot!(u_exact, l=(2, :black), label="reference")

Now we can perform a convergence study.

n = [round(Int, 5 * 10^k) for k in 0:0.5:3]
err = []
for n in n
    t, u = FNC.euler(ivp, n)
    push!(err, norm(u_exact.(t) - u, Inf))
end
@pt :header=["n", "inf-norm error"] [n err]

The error is approximately cut by a factor of 10 for each increase in $n$ by the same factor. A log-log plot also confirms first-order convergence. Keep in mind that since $h=(b-a)/n$ , it follows that $O(h)=O(n^{-1})$ .

plot(n, err;
    m=:o, label="results",
    xaxis=(:log10, L"n"),  yaxis=(:log10, "inf-norm global error"),
    title="Convergence of Euler's method")

# Add line for perfect 1st order.
plot!(n, 0.5 * err[end] * (n / n[end]) .^ (-1), l=:dash, label=L"O(n^{-1})")

Example 6.2.1

We consider the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

As usual, we need to define the function for the right-hand side of the ODE, the interval for the independent variable, and the initial value.

f = @(t, u) sin((t + u)^2);
u0 = -1;
tspan = [0, 4];

Here is the call to Function 6.2.1.

[t, u] = eulersys(f, tspan, u0, 20);
clf, plot(t, u, '.-', displayname='20 steps')
xlabel('t')
ylabel('u(t)')
title(('Solution by Euler''s method'));

[t, u] = eulersys(f, tspan, u0, 50);
hold on, plot(t, u, '.-',  displayname='50 steps')
legend();

Increasing $n$ changed the solution noticeably. Since we know that interpolants and finite differences become more accurate as $h\to 0$ , we should anticipate the same behavior from Euler’s method. We don’t have an exact solution to compare to, so we will use a built-in solver to construct an accurate reference solution.

opt = odeset('RelTol', 1e-13, 'AbsTol', 1e-13);
sol = ode45(f, tspan, u0, opt);
u_ref = @(t) deval(sol, t)';

Now we can perform a convergence study.

n = round(5 * 10.^(0:0.5:3));
err = [];
for k = 1:length(n)
    [t, u] = eulersys(f, tspan, u0, n(k));
    err(k) = norm(u_ref(t) - u, Inf);
end
table(n', err', VariableNames=["n", "inf-norm error"])

clf
loglog(n, err, 'o-')
hold on, loglog(n, 0.5 * err(end) * (n / n(end)).^(-1), '--')
xlabel('n')
ylabel('inf-norm error')
title('Convergence of Euler''s method')
legend('error', 'O(n^{-1})', 'location', 'southwest');

Example 6.2.1

We consider the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

f = lambda t, u: sin((t + u) ** 2)
tspan = [0.0, 4.0]
u0 = -1.0
t, u = FNC.euler(f, tspan, u0, 20)

fig, ax = subplots()
ax.plot(t, u[0, :], "-o", label="$n=20$")
xlabel("$t$"), ylabel("$u(t)$")
title("Solution by Euler's method")
legend();

t, u = FNC.euler(f, tspan, u0, 200)
ax.plot(t, u[0, :], label="$n=200$")
ax.legend()
fig

Increasing $n$ changed the solution noticeably. Since we know that interpolants and finite differences become more accurate as $h\to 0$ , we should anticipate the same behavior from Euler’s method. We don’t have an exact solution to compare to, so we will use solve_ivp to construct an accurate reference solution.

from scipy.integrate import solve_ivp
sol = solve_ivp(f, tspan, [u0], dense_output=True, atol=1e-8, rtol=1e-8)
ax.plot(t, sol.sol(t)[0, :], "--", label="accurate")
ax.legend()
fig

Now we can perform a convergence study.

n_ = array([int(5 * 10**k) for k in arange(0, 3, 0.5)])
err_ = zeros(6)
results = PrettyTable(["n", "error"])
for j, n in enumerate(n_):
    t, u = FNC.euler(f, tspan, u0, n)
    err_[j] = norm(sol.sol(t)[0, :] - u[0, :], inf)
    results.add_row((n, err_[j]))
print(results)

+------+-----------------------+
|  n   |         error         |
+------+-----------------------+
|  5   |   2.7342049884036537  |
|  15  |  0.15019897709239743  |
|  50  |  0.029996197020050186 |
| 158  |  0.008850284724309654 |
| 500  | 0.0027366205261378784 |
| 1581 | 0.0008596857693511373 |
+------+-----------------------+

loglog(n_, err_, "-o", label="results")
plot(n_, 0.5 * (n_ / n_[0])**(-1), "--", label="1st order")
xlabel("$n$"), ylabel("inf-norm error")
title("Convergence of Euler's method")
legend();

Euler’s method is the ancestor of the two major families of IVP methods presented in this chapter. Before we describe them, though, we generalize the initial-value problem itself in a crucial way.

6.2.3Exercises¶

Exercise 6.2.1

✍ Do two steps of Euler’s method for the following problems using the given step size $h$ . Then, compute the error using the given exact solution.

(a) $u' = -2t u, \ u(0) = 2;\ h=0.1;\ \hat{u}(t) = 2e^{-t^2}$

(b) $u' = u + t, \ u(0) = 2;\ h=0.2;\ \hat{u}(t) = -1-t+3e^t$

(c) $t u' + u = 1, \ u(1) = 6, \ h = 0.25;\ \hat{u}(t) = 1+5/t$

(d) $u' - 2u(1-u) = 0, \ u(0) = 1/2, \ h = 0.25; \ \hat{u}(t) = 1/(1 + e^{-2t})$

Exercise 6.2.2

⌨ For each IVP, solve the problem using Function 6.2.1. (i) Plot the solution for $n=320$ . (ii) For $n=10\cdot2^k$ , $k=2,3,\ldots,10$ , compute the error at the final time and make a log-log convergence plot, including a reference line for first-order convergence.

(a) $u' = -2t u, \ 0 \le t \le 2, \ u(0) = 2;\ \hat{u}(t) = 2e^{-t^2}$

(b) $u' = u + t, \ 0 \le t \le 1, \ u(0) = 2;\ \hat{u}(t) = -1-t+3e^t$

(c) $(1+t^3)uu' = t^2,\ 0 \le t \le 3, \ u(0) =1;\ \hat{u}(t) = [1+(2/3)\ln (1+t^3)]^{1/2}$

(d) $u' - 2u(1-u) = 0, \ 0 \le t \le 2, \ u(0) = 1/2; \ \hat{u}(t) = 1/(1 + e^{-2t})$

(e) $v' - (1+x^2) v = 0, \ 1 \le x \le 3, \ v(1) = 1, \ \hat{v}(x) = e^{(x^3+3x-4)/3}$

(f) $v' + (1+x^2) v^2 = 0, \ 0 \le x \le 2, \ v(0) = 2, \ \hat{v}(x) = 6/(2x^3+6x+3)$

(g) $u' = 2(1+t)(1+u^2), \ 0 \le t \le 0.5, \ u(0) = 0, \ \hat{u}(t) = \tan(2t + t^2)$

Footnotes¶

Another point of view is that we can of course make local errors smaller by chopping $h$ in half, but then we have to take twice as many steps. The important quantity, then, is local error per unit step length, which is how $\tau$ is defined.
↩