Accuracy and stability of numerical methods

1.3. Accuracy and stability of numerical methods#

The numerical methods we introduced in the previous section are based on approximations of the derivatives. In this section, we will look closer at the errors introduced by this approximation, how they affect the solution, and how we can estimate their order of magnitude.

Round-off errors#

There are two types of errors: round-off and truncation errors. The round-off error is due to the computer’s limitation to represent a floating number (decimal). To illustrate it consider the following: the difference between 1 and 0.9 is 0.1. If you subtract 0.1 from itself, you should obtain 0, even if this subtraction is repeated ten thousand times. However, in practice, an extremely small error accumulates over each repetition, so the result will not be zero. The following code demonstrates this behavior, using floating-point numbers with varying precision.

def accumulated_error(a, b, solution, iterations):
    error = solution - np.abs(a - b)
    error_accumulated = 0
    for i in range(iterations):
        error_accumulated = error_accumulated + error

    return error_accumulated


print(
    "Error using 16 bits of memory = ",
    accumulated_error(np.float16(1.0), np.float16(0.9), np.float16(0.1), 10000),
)
print(
    "Error using 32 bits of memory = ",
    accumulated_error(np.float32(1.0), np.float32(0.9), np.float32(0.1), 10000),
)
print(
    "Error using 64 bits of memory = ",
    accumulated_error(np.float64(1.0), np.float64(0.9), np.float64(0.1), 10000),
)

As you can see, this simple operation can give discernible errors. Imagine a computation spanning 100 years with a time step of seconds and more complex operations: the round-off error will be present! As you can see, this can be reduced by increasing the precision or the number of digits used to represent numbers but be careful, this is not free as the memory the computer uses increases as well as the computation time.

Truncation errors#

The truncation error is related to the method chosen to approximate the slope. The error can be obtained using our reliable Taylor series expansion which gives the exact solution. Therefore the truncation error per step is:

\[ \text{Truncation error }= y^{TSE}(x_{i+1}) - y^{\text{Forward Euler}}(x_{i+1}) \]

\[ \text{Truncation error }= y(x_i)+\Delta x y'_i +\frac{(\Delta x)^2}{2!} y''_i + ... - y(x_i)-\Delta x y'_i \approx \mathcal{O}(\Delta x^2) \]

The total truncation error accounts for the number of steps as:

\[ \text{Total truncation error } \approx \sum_{i=0}^{n-1} \frac{(\Delta x)^2}{2!} y''_i \approx \frac{(\Delta x)^2}{2!} \frac{b-a}{\Delta x} \bar{y}'' \approx \mathcal{O}(\Delta x) \]

Because this total truncation error is of first order, the forward Euler method is referred to as a first-order method.

Global errors#

The global error is the sum of round-off and total truncation errors. The plot below displays the global errors for both the forward and backward Euler methods, alongside the reference line representing the slope of the truncation error \(\mathcal{O}(\Delta x)\) (1 to 1). As shown, the global errors for both methods align with the expected slope, indicating that the errors decrease in accordance with the methods’ first-order accuracy.

\[ \text{Global Error} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( y_i - y_{\text{exact},i} \right)^2} \]

import numpy as np
import matplotlib.pyplot as plt


def global_error(dx_values):
    errors_forward = []
    errors_backward = []

    for dx in dx_values:
        x = np.arange(0, 30 + dx, dx)
        y_forward = np.zeros(len(x))
        y_backward = np.zeros(len(x))

        alpha = 1
        y_forward[0] = 1
        for i in range(len(x) - 1):
            y_forward[i + 1] = y_forward[i] + dx * (-alpha * y_forward[i])

        y_backward[0] = 1
        for i in range(len(x) - 1):
            y_backward[i + 1] = y_backward[i] / (1 + dx * alpha)

        y_exact = np.exp(-alpha * x)
        n = (30 - 0) / dx
        # Compute the global error (L2 norm) for both methods
        error_forward = np.sqrt((1 / (n - 1)) * np.sum((y_forward - y_exact) ** 2))
        error_backward = np.sqrt((1 / (n - 1)) * np.sum((y_backward - y_exact) ** 2))

        errors_forward.append(error_forward)
        errors_backward.append(error_backward)

    # Plot the global error vs. dx in a loglog plot for both methods
    plt.figure(figsize=(8, 5))
    plt.loglog(
        dx_values, errors_forward, label="Forward Euler global error", marker="o"
    )
    plt.loglog(
        dx_values, errors_backward, label="Backward Euler global error", marker="s"
    )
    plt.loglog(
        dx_values,
        [dx**1 for dx in dx_values],
        label=r"$\mathcal{O} (\Delta x)$",
        linestyle="--",
    )

    plt.title(r"Global error comparison: Forward vs Backward Euler")
    plt.xlabel(r"$\Delta x$")
    plt.ylabel("Global error")
    plt.grid(True)
    plt.legend()
    plt.show()


dx_values = [0.00125, 0.0025, 0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64]
# Call the function to plot the global error for both Forward and Backward Euler
global_error(dx_values)

Stability#

The error that is introduced in each step of the numerical solution ideally does not to increase as the solution advances. Under a well posed and proper solution, the error is expected to reduce with smaller steps. In some cases, the error increases without bound as the solution advances from the initial condition (even with smaller steps): the solution becomes unstable. The stability depends on the numerical method, the step size and the behavior of the differential equation. Therefore, the stability conditions will differ when applying the same numerical method to different equations.

Explicit Euler scheme#

Let’s consider the stability of the forward Euler method applied to a more general form of the problem above.

\[ \frac{dy}{dx} = y'=-\alpha y \]

With initial condition \(y(0)=1\) and \(\alpha>0\), the exact solution is:

\[ y^{\text{exact}}(x_{i})=y_0 e^{-\alpha x_{i}} \]

The forward Euler equivalent is

\[ y_{i+1}=y_{i}-\alpha y_i \Delta x=y_i(1-\alpha \Delta x) \]

Following the initial steps of the numerical solution a pattern arises:

\[ y_{1}=y_0(1-\alpha \Delta x) \]

\[ y_{2}=y_1(1-\alpha \Delta x) =y_0(1-\alpha \Delta x)^2 \]

\[ y_{3}=y_2(1-\alpha \Delta x) =y_0(1-\alpha \Delta x)^3 \]

\[ y_{n}=y_0(1-\alpha \Delta x)^n \]

Comparing this last equation with the exact solution, it can be seen that the term \((1-\alpha \Delta x)^n\) in the numerical solution is approximating the term \(e^{-\alpha \Delta x_i}\) in the exact solution. The latter tends to decay with larger values of \(x_i\) and positive \(\alpha\). To make sure that the term \((1-\alpha \Delta x)^n\) decays with larger \(n\) values \(1-\alpha \Delta x\) should be less than \(|1|\), i.e.,

\[ |1-\alpha \Delta x| < 1 \implies 0 < \alpha \Delta x < 2 \]

Click here for a more detailed derivation

\[ -1 <1-\alpha \Delta x < 1 \]

Subtract 1 of both sides

\[ -2 <-\alpha \Delta x < 0 \]

Divide by -1

\[ 0 <\alpha \Delta x < 2 \]

This is the stability criterion. If we do not comply with it, the solution will be unstable. We can say that the forward Euler scheme is conditionally stable. This is true for every explicit numerical method.

Exercise

Lets go back to the last example described in the code above (for \(\alpha=1\)).

\[ \frac{dy}{dx}= y'=-y \]

We want to know the solution in the domain \([a=0,b=30]\). The initial condition is \(y(x_0)=1\) and the step size is \(\Delta x = 0.2\).

What is the step size at which the problem described using the forward Euler method becomes unstable? Feel free to use the interactive figure below and change \(dx\) to test what happens when the step size becomes larger.

Click rocket –>Live Code to interact with the plot below

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets


def forward_euler(dx):
    # Define the x range and initialize y
    x = np.arange(0, 30 + dx, dx)
    y = np.zeros(len(x))

    # Forward Euler Implementation
    alpha = 1
    y[0] = 1
    for i in range(len(x) - 1):
        y[i + 1] = y[i] + dx * (-alpha * y[i])

    x_exact = np.arange(0, 30 + dx, 0.2)
    # Exact solution
    y_exact = np.exp(-alpha * x_exact)

    # Plot both solutions
    plt.figure(figsize=(8, 5))
    plt.plot(x_exact, y_exact, label="Exact solution", linestyle="--")
    plt.plot(x, y, label="Forward Euler solution", marker="o")
    plt.legend()
    plt.title(f"Forward Euler with dx = {dx}")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid(True)
    plt.show()


# Create an interactive slider for dx
interact(
    forward_euler,
    dx=widgets.FloatSlider(value=0.2, min=0.01, max=3.0, step=0.01, description="dx"),
);

Solution

You can see in the interactive figure above that when the step size becomes larger the forward Euler solution becomes less accurate. You can also see that when the step size becomes bigger than 2 the forward Euler solution starts to blow up. This is because the stability criterion of the forward Euler method is violated.

\[ 0 <\alpha \Delta x < 2 \]

When \(\alpha=1\), then \(\Delta x < 2\) (Note that \(\Delta x\) cannot be negative and therefore will always be greater than 0.)

Note that even when the stability criterion is met, the solution can look “unstable”, for example, show oscillations – as long as it eventually converges. Only when the errors get bigger and bigger over time, we call it unstable.

Implicit Euler scheme#

Exercise:

Modify the step size in the Implicit Euler code, try to make the solution unstable. What do you notice?

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets


def backward_euler(dx):
    # Define the x range and initialize y
    x = np.arange(0, 30 + dx, dx)
    y = np.zeros(len(x))

    ##-----------------------------
    ##Backward Euler Implementation
    ##-----------------------------
    alpha = 1
    y[0] = 1
    for i in range(len(x) - 1):
        y[i + 1] = y[i] / (1 + dx * alpha)
    ##------------------------------

    x_exact = np.arange(0, 30 + dx, 0.2)
    # Exact solution
    y_exact = np.exp(-alpha * x_exact)

    # Plot both solutions
    plt.figure(figsize=(8, 5))
    plt.plot(x_exact, y_exact, label="Exact solution", linestyle="--")
    plt.plot(x, y, label="Backward Euler solution", marker="o")
    plt.legend()
    plt.title(f"Backward Euler with dx = {dx}")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid(True)
    plt.show()

# Create an interactive slider for dx
interact(
    backward_euler,
    dx=widgets.FloatSlider(value=0.2, min=0.01, max=5.0, step=0.01, description="dx"),
);

Let’s understand what happens by following the same procedure as for the stability analysis of explicit Euler.

The exact solution is the same:

\[ y^{\text{exact}}(x_{i})=y_0 e^{-\alpha x_{i}} \]

For the numerical solution, the numerical solution pattern is slightly different:

\[ y_{1}=y_0\left(\frac{1}{1+\alpha \Delta x}\right) \]

\[ y_{2}=y_1\left(\frac{1}{1+\alpha \Delta x}\right) =y_0\left(\frac{1}{1+\alpha \Delta x}\right)^2 \]

\[ y_{3}=y_2\left(\frac{1}{1+\alpha \Delta x}\right) =y_0\left(\frac{1}{1+\alpha \Delta x}\right)^3 \]

\[ y_{n}=y_0\left(\frac{1}{1+\alpha \Delta x}\right)^n \]

This last term to the power \(n\) approximates the decaying exponential in the exact solution, just like in the explicit Euler method. Here, the only condition that must be satisfied for stability is

\[\left|\frac{1}{1+\alpha \Delta x}\right| < 1\]

This can be expressed as:

\[ -1 <\frac{1}{1+\alpha \Delta x} < 1\]

Lets first look at the right side condition \(\frac{1}{1+\alpha \Delta x} < 1\).

If we multiply by the term \(1+\alpha \Delta x\), we get:

\[ 1 < 1+\alpha \Delta x \]

As \(\alpha\) is positive and the step size must be larger than 0, the right side of the criterion is always true.

Lets now look at the left side condition \(-1 <\frac{1}{1+\alpha \Delta x}\):

We again multiply by \(1+\alpha \Delta x\):

\[ -1 - \alpha \Delta x < 1 \]

Again this left side criterion always holds.

Hence, an implicit backward Euler scheme is unconditionally stable. This is also the case for more advanced implicit schemes.

Attribution

This chapter is written by Jaime Arriaga Garcia, Anna Störiko, Justin Pittman and Robert Lanzafame. Find out more here.