Euler’s Method

6.2. Euler’s Method#

Last revised on October 23, 2025, correcting some typos and some reformatting.

References:

import numpy as np
# Shortcuts for some favorite commands:
from numpy import linspace
from matplotlib.pyplot import figure, plot, grid, title, xlabel, ylabel, legend

6.2.1. Test Cases#

A lot of useful intuition comes from solving for the first four examples in the chapter introduction:

  • Example 6.1 \(\displaystyle du/dt = f(t)\) (Integration),

  • Example 6.2 \(\displaystyle du/dt = K u\) (the simplest “genuine” ODE),

  • Example 6.3 \(\displaystyle du/dt = u^2\) (a nonlinear ODE) and

  • Example 6.4 \(\displaystyle du/dt = -\sin t -K(u - \cos t)\) (with fast and slow time scales).

6.2.2. The Tangent Line Method, a.k.a. Euler’s Method#

Once we know \(u(t)\) (or a good approximation) at some time \(t\), we also know the value of \(u'(t) = f(t, u(t))\) there; in particular, we know that \(u(a) = u_0\) and so \(u'(a) = f(a, u_0)\).

This allows us to approximate \(u\) for slightly larger values of the argument (which I will call “time”) using its tangent line:

\[ u(a+h) \approx u(a) + u'(a) h = u_0 + f(a, u_0) h \quad \text{for ``small'' $h$} \]

and more generally

\[ u(t+h) \approx u(t) + f(t, u(t)) h \quad \text{for ``small'' $h$} \]

This leads to the simplest approximation: choose a step size \(h\) determining equally spaced times \(t_i = a + i h\) and define — recursively — a sequence of approximations \(U_i \approx u(t_i)\) with

\[\begin{split}\begin{split} U_0 &= u_0 \\ U_{i+1} &= U_i + h f(t_i, U_i)) \end{split}\end{split}\]

If we choose a number of time steps \(n\) and set \(h = (b-a)/n\) for \(0 \leq i \leq n\), the second equation is needed for \(0 \leq i < n\), ending with \(U_n \approx u(t_n) = u(b)\).

This “two-liner” does not need a pseudo-code description; instead, we can go directly to a rudimentary Python function for Euler’s Method:

def euler_method(f, a, b, u_0, n):
    """Solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0"""
    h = (b-a)/n
    t = linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        u[i+1] = u[i] + f(t[i], u[i])*h
    return (t, u)

See Exercise 6.1.

6.2.3. Examples#

6.2.3.1. Solving for Example 6.1, an integration#

# For integration of -sin(t)

def f_1(t, u):
    """For integration of -sin(x)
    The general solution is
        u(t) = cos(x) + C,
        C = u_0 - cos(a)
    """
    return -np.sin(t)

def u_1(t, a, u_0):
    return np.cos(t) + (u_0 - np.cos(a));

# A helper function for rounding some output to four significant digits
# TO DO: check Python syntax
def approx4(x):
    return round(x, sigdigits=4)

a = 0.0
b = 3/4*np.pi
u_0 = 3.0
n = 20

(t, U) = euler_method(f_1, a, b, u_0, n)  
u = u_1(t, a, u_0)

figure(figsize=[12,6])
title(f"The exact solution is y = cos(x) + {u_0 - np.cos(a)}")
plot(t, u, "g", label="Exact solution")
plot(t, U, ".:b", label=f"Euler's answer for h={(b-a)/n:0.4g}")
legend()
grid(True);
../_images/8aa1264718888995049dc8a96dc12bc06b81ec12530970a2018b7359c9c95301.png

6.2.3.2. Solving for Example 6.2, some exponential functions#

def f_2(t, u):
    """For solving du/dt = K u.
    The variable K may be defined later, so long as that is done before this function is used.
    """
    return K*u

def u_2(t, a, u_0, K):
    return u_0 * np.exp(K*(t-a));

# You could experiment by changing these values here;
# for now I instead redefine them below.
K = 1.0
u_0 = 0.8
a = 0.0
b = 2.0
n = 40

(t, U) = euler_method(f_2, a, b, u_0, n)
u = u_2(t, a, u_0, K)

figure(figsize=[12,6])
title(f"The exact solution is $u = {u_0}\\ \\exp({K}\\ t)$")
plot(t, u, "g", label="Exact solution")
plot(t, U, ".:b", label=f"Euler's answer for h={(b-a)/n:0.4g}")
legend()
grid(True);
../_images/c300043855a4b8845284878c1a7ff2922a4b30b997c33435ba68ca417b4d8012.png 
# You could experiment by changing these values here.
K = -0.5
u_0 = 3.
a = 0.
b = 2.

(t10, u_10) = euler_method(f_2, a, b, u_0, 10)
(t20, U20) = euler_method(f_2, a, b, u_0, 20)
t = t20
u = u_2(t, a, u_0, K)

figure(figsize=[12,6])
title(f"The exact solution is $y = {u_0} \\, \\exp({K} \\, t)$")
plot(t, u, "g", label="Exact solution")
plot(t10, u_10, ".:r", label=f"Euler's answer for h={(b-a)/10:0.4g}")
plot(t20, U20, ".:b", label=f"Euler's answer for h={(b-a)/20:0.4g}")
legend()
grid(True);
../_images/b99c363c39f760d52ec21193239e63da25d4d0190868fab81f4403753a0c3d6e.png

6.2.3.3. Solving for Example 6.3: solutions that blow up#

def f_3(t, u):
    """For solving du/dt = u^2
    # The general solution is u(t) = 1/((a + 1/u_0) - t), = 1/(T-t) with T = a + 1/u_0
    """
    return u**2

def u_3(t, a, u_0):
    T = a + 1.0/u_0
    return 1.0/(T - t);

a = 0.0
b = 0.9
u_0 = 1.0

(t100, U100) = euler_method(f_3, a, b, u_0, 100)
(t200, U200) = euler_method(f_3, a, b, u_0, 200)
t = t200
u = u_3(t, a, u_0)

figure(figsize=[12,6])
title(f"The exact solution is $u = 1/({a + 1/u_0} - t)$")
plot(t, u, "g", label=f"Exact solution")
plot(t100, U100, ".:r", label=f"Euler's answer for h={(b-a)/100:0.2g}")
plot(t200, U200, ".:b", label=f"Euler's answer for h={(b-a)/200:0.2g}")
legend()
grid(True);
../_images/7fa7fc69d140e105ea3ca0fd73cc4de9327d9023193948f960a49b61ecadab85.png

There is clearly a problem when \(t\) reaches 1; let us explore that:

a = 0.0
b = 0.999
u_0 = 1.0
n = 200

(t, U) = euler_method(f_3, a, b, u_0, n)
T = a + 1/u_0
tplot = linspace(a, b, 1000)  # More t values are needed to get a good graph near the vertical asymptote
u = u_3(tplot, a, u_0)

figure(figsize=[12,6])
title(f"The exact solution is $u = 1/({T} - t)$")
plot(tplot, u, "g", label=f"Exact solution")
plot(t, U, ":b", label=f"Euler's answer for h={(b-a)/n:0.4g}")
legend()
grid(True);
../_images/329d6fad325a6bffc91f113d92af888cf749fdc5feb4251f1dc80d5d38bf9bcc.png

Clearly Euler’s method can never produce the vertical asymptote.

The best we can do is improve accuracy by using more, smaller time steps:

n = 10000

(t, U) = euler_method(f_3, a, b, u_0, n)
T = a + 1/u_0
tplot = linspace(a, b, 1000)  # More t values are needed to get a good graph near the vertical asymptote
u = u_3(tplot, a, u_0)

figure(figsize=[12,6])
title(f"The exact solution is $u = 1/({T} - t)$")
plot(tplot, u, "g", label="Exact solution")
plot(t, U, ":b", label=f"Euler's answer for h={(b-a)/n:0.4}")
legend()
grid(True);
../_images/adc967c8af7047a6caab0f596563da25d9873daba7da9fe763c6bac1de891d61.png

6.2.3.4. Solving for Solving for Example 6.4, a stiff ODE#

def f_4(t, u):
    """The variable K may be defined later, so long as that is done before this function is used.
    The general solution is u(t) = u(t; a, u_0, K) = \\cos t + (u_0 - cos(a)) e^{-K (t-a)}
    """
    return -np.sin(t) - K*(u - np.cos(t))

def u_4(t, a, u_0, K):
    return np.cos(t) + (u_0 - np.cos(a)) * np.exp(K*(a-t))

With enough steps (small enough step size \(h\)), all is well:

a = 0.0
b = 2 * np.pi  # One period
u_0 = 2.0
K = 40.0
n = 400

(t, U) = euler_method(f_4, a, b, u_0, n)
u = u_4(t, a, u_0, K)

figure(figsize=[12,6])
title(f"The exact solution is u = cos t + {u_0-1:0.4g} exp(-{K} t)")
plot(t, u, "g", label=f"Exact solution for {K=}")
plot(t, U, ":b", label=f"Euler's answer for h={(b-a)/n:0.4g}")
legend()
grid(True);
../_images/b65b4c3869e76e1cf10780ef6f243af7afcfb04cb83115c53a59714cef8c3430.png

However, with large steps (still small enough to handle the \(\cos t\) part), there is a catastrophic failure, with growing oscillations.

As we will see, these are a characteristic feature of instability.

n = 124

(t, U) = euler_method(f_4, a, b, u_0, n)
u = u_4(t, a, u_0, K)

figure(figsize=[12,6])
title(f"The exact solution is u = cos t + {u_0-1:0.3} exp(-{K} t)")
plot(t, u, "g", label=f"Exact solution for {K=}")
plot(t, U, '.:b', label=f"Euler's answer for h={(b-a)/n:0.3}")
legend()
grid(True);
../_images/33f4e355b3363c69b0c88bf9533b06b9ca518826667724800610627c9cffedf1.png

To show that the \(K\) part is the problem, reduce \(K\) while leaving the rest unchanged:

K = 10.0

(t, U) = euler_method(f_4, a, b, u_0, n)
u = u_4(t, a, u_0, K)

figure(figsize=[12,6])
title(f"The exact solution is u = cos t + {u_0-1:0.3} exp(-{K} t)")
plot(t, u, "g", label=f"Exact solution for {K=}")
plot(t, U, '.:b', label=f"Euler's answer for h={(b-a)/n:0.3}")
legend()
grid(True);
../_images/0717fc3a509d0a3913cdedff1efe48031a9dde610685fe29adb74b115c754549.png

6.2.4. Variable Time Step Sizes \(h_i\) (just a preview)#

It is sometime useful to adjust the time step size; for example reducing it when the derivative is larger, (as happens in Example 3 above). This gives a slight variant, now expressed in pseudo-code:

Input: \(f\), \(a\), \(b\), \(n\)
\(t_0 = a\)
\(U_0 = u_0\)
for i from 1 to n
\(\quad\) Choose \(h_i\) somehow
\(\quad t_{i} = t_{i-1} + h_i\)
\(\quad U_{i} = U_{i-1} + h_i f(t_{i-1}, U_{i-1})\)
end

In section Error Control and Variable Step Sizes we will see how to estimate errors within an algorithm, and then how to use such error estimates to guide the choice of step size.

6.2.5. Error Analysis for the Canonical Test Case, \(u' = K u\).#

A great amount of intuition about numerical methods for solving ODE IVPs comes from that “simplest nontrivial example”, number 2 above. We can solve it with constant step size \(h\), and thus study its errors and accuracy. The recursion relation is now

\[U_{i} = U_{i-1} + h K U_{i-1} = U_{i-1} (1 + h K),\]

with solution

\[U_i = u_0 (1 + K h)^i\]

For comparison, the exact solution of this ODE IVP is

\[u(t_i) = u_0 e^{K(t_i - a)} = u_0 e^{K ih} = u_0 (e^{K h})^i\]

So each is a geometric series: the difference is that the growth factor is \(G = (1 + K h)\) for Euler’s method, vs \(g = e^{K h} = 1 + K h + (K h)^2/2 + \cdots = 1 + K h + O(h^2)\) for the ODE.

Ths deviation at each time step is \(O(h^2)\), suggesting that by the end \(t=b\), at step \(n\), the error will be \(O(n h^2) = O\left(\displaystyle\frac{b-a}{h} h^2\right) = O(h)\).

This is in fact what happens, but to verify that, we must deal with the challenge that once an error enters at one step, it is potentially amplified at each subsequent step, so the errors introduced at each step do not simply get summed like they did with definite integrals.

6.2.5.1. Global Error and Local (Truncation) Error#

Ultimately, the error we need to understand is the global error: at step \(i\),

\[ E_i = u(t_i) - U_i \]

We will approach this by first considering the new error added at each step, the local truncation error (or discretization error).

At the first step this is the same as above:

\[ e_1 = u(t_1) - U_1 = u(a+h) - U_1 \]

However at later steps we compare the results \(U_{i+1}\) to what the solution would be if it were exact at the start of that step: that is, if \(U_i\) were exact.

Using the notation \(u(t; t_i, U_i)\) introduced above for the solution of the ODE with initial condition \(u(t_i) = U_i\), the location truncation error at step \(i\) is the discrepancy at time \(t_{i+1}\) between what Euler’s method and the exact solution give when both start at that point \((t_i, U_i)\):

\[ e_i = u(t; t_i, U_i) - U_{i+1} \]

6.2.5.2. Error propagation in \(u' = K u\), \(K \geq 0\).#

After one step, \(E_1 = u(t_1) - U_1 = e_1\).

At step 2,

\[ E_2 = u(t_2) - U_2 = (u(t_2) - u(t_2, t_1, U_1) + (u(t_2, t_1, U_1) - U_2) = (u(t_2) - u(t_2, t_1, U_1) + e_2 \]

The first term is the difference at \(t=t_2\) of two solutions with values at \(t = t_1\) being \(u(t_1)\) and \(U_1\) respectively. As the ODE is linear and homogeneous, this is the solution of the same ODE with value at \(t=t_1\) being \(u(t_1) - U_1\), which is \(e_1\): that solution is \(e_1 e^{K(t - t_1)}\), so at \(t = t_2\) it is \(e_1 e^{K h}\). Thus the global error after two steps is

\[ E_2 = e_2 + (e^{K h})e_1: \]

the error from the previous step has been amplified by the growth factor \(g = e^{K h}\):

\[ E_2 = e_2 + g e_1: \]

This continues, so that

\[ E_3 = e_3 + g E_2 = e_3 + g (e_2 + g e_1) = e_3 + g e_2 + g^2 e_1 \]

and so on, leading to

\[ E_i = e_i + g e_{i-1} + g^2 e_{i-2} + \cdots + g^{i-1} e_{1}. \]

6.2.5.3. Bounding the local truncation errors …#

To get a bound on the global error from the formula above, we first need a bound on the local truncation errors \(e_i\).

Taylor’s theorem gives \(e^{K h} = 1 + K h + e^{K \xi} (K h)^2/2\), \(0 \leq \xi \leq h\), so

\[ e_i = U_i e^{K h} - U_i (1 + K h) = U_i (e^{K \xi} h^2/2) \]

Since \(0 \leq \xi \leq h\) we have \(1 \leq e^{K \xi} \leq e^{K h}\), giving

\[|e_i| \leq |U_i| \frac{e^{K h}}{2} h^2\]

Also, since \(1 + K h < e^{K h}\), \(|U_i| < |u(t_i)| = |u_0| e^{K(t_i - a)},\) and we only need this to the beginning of the last step, \(i \leq n-1\), for which

\[ |U_i| < |u_0| e^{K(b - h - a)} \]

Thus

\[ |e_i| \leq \frac{|u_0| e^{K(b - h - a)} e^{K h}}{2} h^2 = \frac{|u_0 e^{K(b - a)}|}{2} h^2 \]

That is,

\[ |e_i| \leq C h^2 \quad \text{where}\; C := \frac{|u_0 e^{K(b - a)}|}{2} \]

6.2.5.4. … and using this to complete the bound on the global truncation error#

Using this bound on the local errors \(e_i\) in the above sum for the global error \(E_i\),

\[ | E_i | \leq C h^2 (1 + g + \cdots g^{i-1}) = C \frac{g^i - 1}{g-1} h^2 \]

Since \(g^i = e^{K h i} = e^{K (t_i - a)}\) and the denominator \(g - 1 = e^{K h} - 1 > K h\), we get

\[ | E_i | \leq C \frac{e^{K (t_i - a)} - 1}{K h} h^2 \leq \frac{|u_0 e^{K(b - a)}|}{2} \frac{e^{K (t_i - a)} - 1}{K} h, = O(h) \]

Note that this global error formula is built from three factors:

  • The first is the constant \(\displaystyle \frac{|u_0 e^{K(b - a)}|}{2}\) which is roughly half of the maximum value of the exact solution over the interval \([a, b]\).

  • The second \(\displaystyle \frac{e^{K (t_i - a)} - 1}{K}\) depends on \(t\) but is bounded by \(\displaystyle \frac{e^{K (b - a)} - 1}{K}\) on the domain \([a,b]\), and

  • The third is \(h\), showing the overall order of accuracy: the overall absolute error is \(O(h)\), so first order.

6.2.5.5. A more general error bound#

A very similar result applies to the solution \(u(t; a, u_0)\) of the more general initial value problem

\[ \frac{du}{dt} = f(t, u), \quad u(a) = u_0 \]

so long as the function \(f\) is “somewhat well-behaved” in that it satisfies a so-called Lipschitz Condition: that there is some constant \(K\) such that

\[ \left| \frac{\partial f}{\partial u}(t, u) \right| \leq K \]

for the relevant time values \(a \leq t \leq b\).

(Aside: As you might have seen in a course on differential equations, such a Lipschitz condition is necessary to even guarantee that the initial value problem has a unique solution, so it is a quite reasonable requirement.)

Then this constant \(K\) plays the part of the exponential growth factor \(K\) above:

first one shows that the local trunction error is bounded by

\[ |e_i| \leq C h^2 \text{ where now } C := \frac{|u_0 e^{K(b - a)}|}{2}; \]

then calculating as above bounds the global truncation error with

\[| E_i | \leq \frac{|u_0 e^{K(b - a)}|}{2} \frac{e^{K (t_i - a)} - 1}{K} h, = O(h)\]

Remark 6.3

We will see a similar result for the large class of so-call one-step methods: if \(|e_i| \leq C h^{p+1}\) then $\(| E_i | \leq C \frac{e^{K (t_i - a)} - 1}{K} h^p\)$ so that the hard part is getting suitable control of the local truncation error.

Remark 6.4

As with definite integrals, first order accuracy is not very impressive, so in the next section on Runge-Kutta Methods we will explore several widely-used methods that improve to second order and then fourth order accuracy. Later, we will see how to get even higher orders.

6.2.6. Error Propagation, for Integration vs a True ODE#

We can illustrate how this exponential “compunding” of errors looks in some examples, and compare to the better-behaved errors in definite integrals.

This will be done by looking at the effect of a small change in the initial value, to simulate an error that arises there.

6.2.6.1. Error propagation for integration: Example 6.1#

a = 0.0
b = 2*np.pi
u_0 = 1.0  # Original value
n = 100

(t, U) = euler_method(f_1, a, b, u_0, n)

But now “perturb” the initial value in all cases by this much:

delta_u_0 = 0.1

(x, U_perturbed) = euler_method(f_1, a, b, u_0+delta_u_0, n)

u = u_1(t, a, u_0);

figure(figsize=[12,6])
title("The solution before perturbing $u(0)$ was $u = \\cos(x)$")
plot(t, u, "g", label="Original exact solution")
plot(t, U, ".:b", label="Euler's answer before perturbation")
plot(t, U_perturbed, ".:r", label="Euler's answer after perturbation")
legend()
grid(True);
../_images/bbb49d17e6775de58af92579ae7496299e48189fbdf88d19d354b5b29658ab71.png

This just shifts all the \(u\) values up by the perturbation of \(u_0\).

6.2.6.2. Error propagation for a “genuine” ODE: Example 6.2#

K = 1.0
a = 0.0
b = 2.0
u_0 = 1.0  # Original value
delta_u_0 = 0.1
n = 100

(t, U) = euler_method(f_2, a, b, u_0, n)
(t, U_perturbed) = euler_method(f_2, a, b, u_0+delta_u_0, n)
u = u_2(t, a, u_0, K)

figure(figsize=[12,6])
title(r"The solution before perturbing $u(0)$ was $u = $"+f"{u_0}"+r"$\,\exp($"+f"{K}"+r"$ \, t)$")
plot(t, u, "g", label="Original exact solution")
plot(t, U, ".:b", label="Euler's answer before perturbation")
plot(t, U_perturbed, ".:r", label="Euler's answer after perturbation")
legend()
grid(True);
../_images/53247f2ddf372bdc5eca22d48b2491e336b772354e78bfdb10b1902e011197c4.png

Graphing the error shows its exponential growth:

figure(figsize=[12,6])
title("Error")
plot(t, u - U_perturbed, '.:')
U_perturbed
grid(True);
../_images/0cdfc894508682308fa70eb24cdd9d9db9c6e6793d54d5d4d3c6756a6a74c4f5.png

6.2.7. Exercises#

Exercise 6.1

Show that for the integration case \(f(t, u) = f(t)\), Euler’s Method is the same as the composite left-hand endpoint rule, as in the section Definite Integrals, Part 2.

Contents