Finite differences: from equations to computation

Quick refresher

A difference replaces a derivative with an error that the Taylor series shows. The forward and backward differences are of first order: half the step gives half the error. The central difference and the three-point second difference are of second order: the error becomes four times smaller.
The truncation error is the residual (the amount by which the equation fails to hold) that the exact solution gives when it is substituted into the difference equation. A scheme approximates the equation (it is consistent with it) if the residual tends to zero as the grid and the step are refined.
Stability is the property of not amplifying errors from step to step. Approximation together with stability gives convergence, and the error of the solution falls with the same order as the truncation error.
The von Neumann analysis splits the error into harmonics (waves of different lengths) and finds the amplification factor $G$ of each one over a step (the number by which that wave is multiplied in one step). A scheme is stable if $|G| \le 1$ for all harmonics. For the explicit scheme of the pressure equation this gives $r = \eta\Delta t/\Delta x^2 \le 1/2$ .
The implicit scheme (first order in time) is stable at any step and does not oscillate. The Crank–Nicolson scheme is also stable and is more accurate (second order). But with a large step its fast components hardly decay and change sign at every step: after abrupt changes the solution “rings” (it keeps swinging up and down).
A convergence study separates the grid error and the step error. Extrapolation uses the results on several grids to predict the answer for an infinitely fine grid. At a fixed step it gives the answer for that step and not the exact solution, so we check the step separately.

Key relationships

\frac{p_{i+1} - p_{i-1}}{2\Delta x} = p'_i + \frac{\Delta x^2}{6}\,p'''_i + O(\Delta x^4)

central difference: second order

\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2} = p''_i + \frac{\Delta x^2}{12}\,p^{(4)}_i + O(\Delta x^4)

three-point second difference: second order

\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\left[\sigma\,\delta^2 p_i^{n+1} + (1 - \sigma)\,\delta^2 p_i^{n}\right]

weighted scheme; $\delta^2 p_i = (p_{i+1} - 2p_i + p_{i-1})/\Delta x^2$ ; $\sigma$ = 0 is explicit, 1 is implicit, ½ is Crank–Nicolson

G(\xi) = \frac{1 - 4(1 - \sigma)\,r\sin^2(\xi/2)}{1 + 4\sigma r\sin^2(\xi/2)},\qquad r = \frac{\eta\,\Delta t}{\Delta x^2}

amplification factor of a harmonic with phase $\xi$ per cell; stability means $|G| \le 1$

e \approx C_x\,\Delta x^{p} + C_t\,\Delta t^{q}

error model for a convergence study: $p$ and $q$ are the orders in space and in time

Common pitfalls

Taking the order of a scheme for a guarantee of accuracy: it says how fast the error falls when we refine the grid, not how large it is on the grid we chose.
Confusing the convergence of a scheme (the solution tends to the exact one as we refine) with the convergence of the iterations within a step (the equations of the step are solved to a given residual).
Judging stability from the first steps: an unstable component grows out of round-off errors (tiny errors from rounding numbers in the computer) and becomes visible only after tens of steps.
Using the Crank–Nicolson scheme with a large step right after a well is started or shut in: the pressure at the well will oscillate.
Checking convergence only in the grid at a fixed time step, or estimating the order from two calculations instead of three.

After this lesson you can

derive difference formulas for the first and second derivatives from the Taylor series, and find their truncation error and order
write the explicit scheme (flows from the pressures at the start of the step), the implicit scheme (flows from the pressures at the end) and the Crank–Nicolson scheme (flows from the average of the two) for the pressure equation as special cases of a weighted scheme (flows from a weighted mix of the two), and give the order of each in time and in space
tell approximation, stability and convergence apart, and explain why convergence needs the first two together
carry out a von Neumann stability analysis: find the amplification factor of a harmonic (the number by which a wave-shaped error is multiplied in one step), and from it the largest step of the explicit scheme and the properties of the implicit and Crank–Nicolson schemes
plan a convergence study in the grid and in the time step, estimate the observed order, and separate the grid error from the time-step error

From derivatives to differences

In a simulation model the pressure is known not everywhere but only at separate points, the cell centres, and only at separate moments of time. We say it is known on a grid in space and in time. Let $p_i^n$ be the pressure at the point $x_i$ at the moment $t^n$ . Neighbouring points are a step $\Delta x$ apart, and neighbouring moments are a step $\Delta t$ apart. A grid has no derivatives, only numbers at the nodes (the points of the grid). The finite-difference method replaces each derivative by a combination of these numbers, a finite difference. The set of nodes that enter a difference formula is called the stencil.

In the lesson Pressure equation for single-phase flow the equations of the cells were obtained in a different way, from the mass balance of a cell. For equal cells of a uniform formation both ways lead to the same result: the link coefficient divided by the storage capacity of the cell equals the hydraulic diffusivity divided by the square of the cell size. So everything the theory of finite-difference schemes says about error and stability also applies to the cell equations of a simulator.

The cell equation as a finite-difference scheme

S\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = a\left(p_{i+1} - 2p_i + p_{i-1}\right)\quad\Longleftrightarrow\quad \frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\,\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2}

where:

$S = V_p c_t$ storage capacity of the cell, $V_p = A\,\Delta x\,\phi$ , m³/MPa
$a = 0.0853\,kA/(\mu\Delta x)$ link coefficient of two equal neighbouring cells, m³/(d·MPa)
$\eta = 0.0853\,k/(\phi\mu c_t)$ hydraulic diffusivity (how fast pressure changes spread through the formation); $\eta/\Delta x^2 = a/S$ , m²/d
$p_i^{n},\ p_i^{n+1}$ pressure in cell $i$ at the start and at the end of the step, MPa
$\Delta x,\ \Delta t,\ A$ cell size, time step, face area, m, d, m²

Units: metric oilfield units: $k$ in mD, $\mu$ in mPa·s, $c_t$ in 1/MPa, $p$ in MPa, $t$ in d; 0.0853 is the constant of Darcy's law in these units (Rock permeability)

Holds when

Equal cells of a uniform formation, the linear pressure equation, an inner cell without a well.
The moment whose pressures give the right-hand side has not been chosen yet. The scheme depends on it (see below).

On the right is the central second difference multiplied by the diffusivity: the cell equation is a finite-difference scheme for the diffusivity equation. For a heterogeneous formation and cells of different sizes the coefficients change from cell to cell, but the structure is the same.

The Taylor series: truncation error and order

The Taylor series shows how far a difference is from a derivative. It expresses the value of a smooth function at a neighbouring point through the value and the derivatives at the given point.

The Taylor series for neighbouring nodes

p_{i\pm1} = p_i \pm \Delta x\,p'_i + \frac{\Delta x^2}{2}\,p''_i \pm \frac{\Delta x^3}{6}\,p'''_i + \frac{\Delta x^4}{24}\,p^{(4)}_i \pm \ldots

where:

$p_i,\ p_{i\pm1}$ values of the function at node $i$ and at the neighbouring nodes, MPa
$p'_i,\ p''_i,\ p'''_i,\ p^{(4)}_i$ first to fourth derivatives with respect to $x$ at node $i$ , MPa/m, MPa/m², …
$\Delta x$ grid step, m

Units: any consistent units

Holds when

The function is smooth: it has continuous derivatives of the needed order between the nodes.
The grid step is constant (a non-uniform grid is in the advanced section).

The upper signs are for node $i + 1$ , and the lower signs are for $i - 1$ . The smaller the step, the faster the terms of the series fall, and the smaller the part that we discard.

If we subtract $p_i$ and divide by the step, we get the forward difference. If we subtract the two expansions from each other, we get the central difference. If we add them, we get the second difference. The terms that remain after the derivative are the error of the formula.

Difference formulas and their error

\begin{aligned} &\frac{p_{i+1} - p_i}{\Delta x} = p'_i + \frac{\Delta x}{2}\,p''_i + O(\Delta x^2) &&\text{(forward)}\\ &\frac{p_i - p_{i-1}}{\Delta x} = p'_i - \frac{\Delta x}{2}\,p''_i + O(\Delta x^2) &&\text{(backward)}\\ &\frac{p_{i+1} - p_{i-1}}{2\Delta x} = p'_i + \frac{\Delta x^2}{6}\,p'''_i + O(\Delta x^4) &&\text{(central)}\\ &\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2} = p''_i + \frac{\Delta x^2}{12}\,p^{(4)}_i + O(\Delta x^4) &&\text{(second difference)} \end{aligned}

where:

$O(\Delta x^m)$ terms of order $\Delta x^m$ : when the step is halved, they decrease by about $2^m$ times, —
$p'_i,\ p''_i,\ p'''_i,\ p^{(4)}_i$ derivatives at node $i$ , MPa/m, MPa/m², …

Units: any consistent units

Holds when

A smooth function and a constant step.
The values at the nodes are exact (no noise or round-off).

When we subtract the expansions, the terms with even powers of the step cancel, and when we add them, the terms with odd powers cancel. So the error of the central difference and of the second difference starts with the second power of the step, while for the forward and backward differences it starts with the first power. The forward and backward differences are wrong by an amount proportional to the step and to the curvature of the function.

Truncation error (approximation error): The difference between a difference formula and the derivative that it replaces, when the exact function is substituted (truncation error). Here we denote it by $\psi$ .
Order of accuracy: The exponent of the step in the leading term of the error: a formula has order $m$ if $\psi = O(\Delta x^m)$ . When the step is halved, the error falls about $2^m$ times.
Leading term: The first non-vanishing term of the series in the error, for example $\tfrac{\Delta x}{2}p''$ for the forward difference. It shows not only the order but also the sign of the error: for $p'' < 0$ the forward difference underestimates the derivative.

Difference formulas for derivatives

Formula: forward, $p'$
Nodes: $i,\ i + 1$
Leading term of the error: $\tfrac{\Delta x}{2}\,p''$
Order: 1

Formula: backward, $p'$
Nodes: $i - 1,\ i$
Leading term of the error: $-\tfrac{\Delta x}{2}\,p''$
Order: 1

Formula: central, $p'$
Nodes: $i \pm 1$
Leading term of the error: $\tfrac{\Delta x^2}{6}\,p'''$
Order: 2

Formula: second difference, $p''$
Nodes: $i,\ i \pm 1$
Leading term of the error: $\tfrac{\Delta x^2}{12}\,p^{(4)}$
Order: 2

Formula	Nodes	Leading term of the error	Order
forward, $p'$	$i,\ i + 1$	$\tfrac{\Delta x}{2}\,p''$	1
backward, $p'$	$i - 1,\ i$	$-\tfrac{\Delta x}{2}\,p''$	1
central, $p'$	$i \pm 1$	$\tfrac{\Delta x^2}{6}\,p'''$	2
second difference, $p''$	$i,\ i \pm 1$	$\tfrac{\Delta x^2}{12}\,p^{(4)}$	2

A uniform grid and a smooth function. Formulas in time are built in the same way: the difference $(p^{n+1} - p^n)/\Delta t$ is a forward difference if we attribute it to the moment $t^n$ , and a central difference if we attribute it to the middle of the step.

Worked examplePressure derivatives near a well from three points

Given (assumed data)

Pressure profile (the well of the lesson Flow regimes): $p(r) - p_{wf} = m\ln(r/r_w)$ ; $m$ = 0.388 MPa = 388 kPa; $r_w$ = 0.1 m
Point and steps: $r$ = 10 m; $\Delta r$ = 2 m and 1 m
Exact derivatives: $p' = m/r$ = 38.8 kPa/m; $p'' = -m/r^2$ = −3.88 kPa/m²

Solution

Pressure differences for a step of 2 m: $p(12) - p(10) = 388\ln 1.2 = 70.741$ kPa; $p(10) - p(8) = 388\ln 1.25 = 86.580$ kPa.
Forward difference: $70.741/2 = 35.37$ kPa/m, error $-3.43$ kPa/m (−8.8%). Backward: $86.580/2 = 43.29$ kPa/m, error $+4.49$ kPa/m. Central: $(70.741 + 86.580)/4 = 39.33$ kPa/m, error $+0.530$ kPa/m (+1.37%).
Step of 1 m: $p(11) - p(10) = 36.980$ kPa, $p(10) - p(9) = 40.880$ kPa. Forward: 36.98 kPa/m (error $-1.82$ ), backward: 40.88 (error $+2.08$ ), central: 38.93 kPa/m (error $+0.130$ , 0.34%).
How many times the errors fell when the step was halved: for the forward difference $3.43/1.82 = 1.9$ , first order; for the central difference $0.530/0.130 = 4.1$ , second order.
Leading terms of the Taylor series for a step of 2 m: for the forward difference $\tfrac{\Delta r}{2}p'' = -3.88$ kPa/m, for the central difference $\tfrac{\Delta r^2}{6}p''' = +0.517$ kPa/m (here $p''' = 2m/r^3$ ). They are close to the actual errors. The remainder comes from the next terms of the series.
Second derivative from three points: for a step of 2 m $(70.741 - 86.580)/2^2 = -3.960$ kPa/m² (error 2.1%), for a step of 1 m $(36.980 - 40.880)/1^2 = -3.900$ kPa/m² (0.5%). The error also fell about four times.

Answer: With a step of 1 m the forward difference is wrong by 4.7%, and the central one by 0.34%. When the step was halved, the error of the first fell 1.9 times, and that of the second fell 4.1 times.

What it means. The central difference with a step of 2 m is more accurate than the forward difference with a step of 1 m: the gain from the higher order is larger than the gain from the smaller step. The leading term predicts the sign of the errors. The curvature of the profile is negative, so the forward difference underestimates the slope and the backward difference overestimates it. Near a well the pressure derivatives are large and grow fast towards the wellbore (the third derivative goes as $1/r^3$ ). So a uniform grid is inefficient here. We make the grid finer near a well (Grids in reservoir simulation models), and we move the pressure drop at the wellbore itself into the well model (Well model in a grid block). The data in this example are assumed for illustration.

Schemes for the pressure equation and their error

Take the one-dimensional linear pressure equation for a uniform formation, the diffusivity equation $\partial p/\partial t = \eta\,\partial^2 p/\partial x^2$ (Diffusivity equation). We replace the second derivative with respect to $x$ by a second difference. What remains is to decide the pressures of which moment to put into it. The lesson Pressure equation for single-phase flow covered the two extreme choices: the start of the step (the explicit scheme) and the end of the step (the implicit scheme). It is convenient to combine them in one weighted scheme: the right-hand side is a weighted average of the second differences at the start and at the end of the step.

The weighted scheme for the pressure equation

\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\left[\sigma\,\frac{p_{i+1}^{n+1} - 2p_i^{n+1} + p_{i-1}^{n+1}}{\Delta x^2} + (1 - \sigma)\,\frac{p_{i+1}^{n} - 2p_i^{n} + p_{i-1}^{n}}{\Delta x^2}\right]

where:

$\sigma$ weight of the end of the step: 0 is the explicit scheme, 1 is the implicit scheme, ½ is the Crank–Nicolson scheme, —
$\eta$ hydraulic diffusivity, $0.0853\,k/(\phi\mu c_t)$ , m²/d
$p_i^{n},\ p_i^{n+1}$ pressure at node (cell) $i$ at the start and at the end of the step, MPa
$\Delta x,\ \Delta t$ grid step and time step, m, d

Units: metric oilfield units (any consistent units)

Holds when

A linear equation, a uniform formation, a uniform grid. For a heterogeneous formation the same is written for the cell equations with the coefficients $a_{ij}$ and $S_i$ .
For $\sigma > 0$ each equation has three unknowns at the end of the step, and at each step we solve a tridiagonal system, as in the implicit scheme.

The weight $\sigma$ sets the moment whose pressures give the flows during the step. For $\sigma = \tfrac12$ the flows are calculated from the average pressure of the step. This is the intermediate scheme that the lesson on the pressure equation mentions.

We find the truncation error of a scheme in the same way as for single formulas. We substitute the exact solution into it and expand all the values in a Taylor series around the middle of the step (the derivation is in the advanced section). The result is:

Truncation error of the weighted scheme

\psi_i^{n} = \left(\tfrac12 - \sigma\right)\Delta t\,\frac{\partial^2 p}{\partial t^2} - \frac{\eta\,\Delta x^2}{12}\,\frac{\partial^4 p}{\partial x^4} + O\!\left(\Delta t^2 + \Delta x^4\right)

where:

$\psi_i^{n}$ truncation error: the left side of the scheme minus the right side when the exact solution is substituted, MPa/d
$\partial^2 p/\partial t^2,\ \partial^4 p/\partial x^4$ derivatives of the exact solution in the middle of the step, MPa/d², MPa/m⁴
$\sigma,\ \eta,\ \Delta t,\ \Delta x$ weight, diffusivity, steps, —, m²/d, d, m

Units: any consistent units

Holds when

A smooth solution (far from sudden changes in how the wells are run, and from the wells themselves).
Constant coefficients, a uniform grid.

The term with $\Delta t$ vanishes only for $\sigma = \tfrac12$ . So the explicit and implicit schemes are of first order in time and second order in space, and the Crank–Nicolson scheme is of second order in both. The error in space is the same for all three: it comes from the second difference.

Three schemes as special cases of the weighted scheme

Scheme: explicit
$\sigma$: 0
Order in $t$: 1
Order in $x$: 2
Stability: for $r \le 1/2$

Scheme: implicit
$\sigma$: 1
Order in $t$: 1
Order in $x$: 2
Stability: for any $r$

Scheme: Crank–Nicolson
$\sigma$: ½
Order in $t$: 2
Order in $x$: 2
Stability: for any $r$

Scheme	$\sigma$	Order in $t$	Order in $x$	Stability
explicit	0	1	2	for $r \le 1/2$
implicit	1	1	2	for any $r$
Crank–Nicolson	½	2	2	for any $r$

$r = \eta\Delta t/\Delta x^2$ is the mesh Fourier number (the step measured in units of the time that the pressure needs to level out over one cell). The stability conditions are derived below by the von Neumann analysis. For a non-uniform grid, for the explicit scheme it is more convenient to use the condition with the characteristic times of the cells from the lesson Pressure equation for single-phase flow.

Approximation, stability, convergence

Approximation (consistency): The truncation error $\psi$ tends to zero when the steps $\Delta x$ and $\Delta t$ tend to zero. The scheme approximates exactly the equation that we need to solve.
Stability: Errors introduced at any step (in the initial data, by round-off, by an inexact solution of the equations of the step) do not grow without limit in the further calculation. Over the finite time of the calculation they stay smaller than their initial size multiplied by a constant number, however small the steps are (numerical stability).
Convergence: As the grid and the step are refined, the solution of the scheme tends to the exact solution of the equation at every point and at every moment.

The link between them is visible from the equation for the error of the solution, $e = p_{\text{exact}} - p$ . The exact solution satisfies the scheme with a residual $\psi$ , and the calculated one satisfies it without a residual. If we subtract one from the other, we get the same scheme for the error, but with a “source” $\psi$ at every step:

Accumulation of error in a stable scheme

\begin{aligned} &e^{n+1} = \mathbf{A}\,e^{n} + \Delta t\,\mathbf{B}\,\psi^{n},\qquad \mathbf{B} = \left(\mathbf{I} - \sigma\Delta t\,\mathbf{L}\right)^{-1},\qquad \mathbf{A} = \mathbf{B}\left(\mathbf{I} + (1 - \sigma)\Delta t\,\mathbf{L}\right)\\ &\left\|e^{n}\right\| \le \left\|e^{0}\right\| + t^{n}\max_k\left\|\psi^{k}\right\| \end{aligned}

where:

$e^{n},\ \psi^{n}$ vectors of the solution error and of the truncation error over all nodes at step $n$ , MPa, MPa/d
$\mathbf{A}$ operator of the step: what the scheme turns the error at the start of the step into, —
$\mathbf{L}$ matrix of second differences multiplied by the diffusivity, 1/d
$\mathbf{I}$ identity matrix, —
$\|\cdot\|,\ t^{n}$ norm of a vector (a measure of its size, for example the root-mean-square) and the time from the start of the calculation, —, d

Units: any consistent units

Holds when

A linear problem. The estimate on the right holds if $\|\mathbf{A}\| \le 1$ (the scheme is stable) and $\|\mathbf{B}\| \le 1$ , which is the case for the weighted scheme for the pressure equation.
For an unstable scheme $\|\mathbf{A}\| > 1$ , and the contribution of each step is multiplied by $\|\mathbf{A}\|$ at all the following steps.

In a stable scheme the error after the time $t^n$ is not larger than the initial error plus the time multiplied by the largest residual. The residual falls as we refine (approximation), so the error falls too (convergence), and with the same order.

This estimate illustrates the main theorem of the theory of finite-difference schemes. For a linear problem that is well posed (it has one solution that depends smoothly on the data), and for a scheme that approximates it, stability is necessary and sufficient for convergence (it is needed, and it is enough). This is the Lax equivalence theorem. Both halves matter. An unstable scheme does not converge even if it approximates the equation: round-off errors of the order of $10^{-15}$ grow by a fixed factor at every step, and sooner or later they fill the whole result. A stable scheme without approximation may converge, but to the solution of a different equation (an example is in the advanced section).

Von Neumann stability analysis

How can we check that $\|\mathbf{A}\| \le 1$ without computing the matrices? For a linear equation with constant coefficients on a uniform grid we can split the error into harmonics, which are sine waves of different wavelengths, and follow each one separately. The scheme turns a harmonic into the same harmonic multiplied by a number $G$ , the amplification factor. John von Neumann proposed the method in the 1940s. It is also called the method of harmonics or the spectral criterion of stability, a test that works on the separate wave components of the error (von Neumann stability analysis).

Let the error at node number $j$ be $G^n\cos(j\xi)$ , where $\xi$ is the phase shift of the harmonic from one node to the next. For small $\xi$ the harmonic is smooth (a long wave). For $\xi = \pi$ it is a “checkerboard”: neighbouring values are equal in size and opposite in sign, and a shorter wave cannot be drawn on the grid. The second difference of such a harmonic is the same harmonic multiplied by $-4\sin^2(\xi/2)/\Delta x^2$ (this follows from the formula for the sum of cosines). If we put the harmonic into the weighted scheme, we get:

Amplification factor of the weighted scheme

G(\xi) = \frac{1 - 4(1 - \sigma)\,r\sin^2(\xi/2)}{1 + 4\sigma\,r\sin^2(\xi/2)},\qquad r = \frac{\eta\,\Delta t}{\Delta x^2}

where:

$G(\xi)$ how many times the amplitude of the harmonic changes over a step; $G < 0$ means that the harmonic changes sign at every step, —
$\xi$ phase shift of the harmonic per cell, $0 < \xi \le \pi$ ; the wavelength is $2\pi\Delta x/\xi$ , rad
$r$ mesh Fourier number: the time step in units of the time that the pressure needs to level out over one cell, —
$\sigma$ weight of the scheme, —

Units: dimensionless

Holds when

A linear equation with constant coefficients and a uniform grid. Boundaries are not taken into account (an infinite or periodic grid).
For a heterogeneous formation we check the condition locally, with the local $\eta$ and $\Delta x$ . For a nonlinear equation we check it with the coefficients “frozen” at the step. The condition is then necessary but not always sufficient.

The exact solution reduces the same harmonic over a step by a factor of $\e^{-r\xi^2}$ . A scheme is stable if $|G| \le 1$ for all $\xi$ , and it is accurate for the harmonics for which $G$ is close to $\e^{-r\xi^2}$ .

Stability conditions of the weighted scheme

\sigma = 0:\ \ r \le \tfrac12;\qquad 0 < \sigma < \tfrac12:\ \ r \le \frac{1}{2(1 - 2\sigma)};\qquad \sigma \ge \tfrac12:\ \ \text{for any } r

where:

$r$ mesh Fourier number $\eta\Delta t/\Delta x^2$ , —
$\sigma$ weight of the scheme, —

Units: dimensionless

Holds when

They follow from $G \ge -1$ for the shortest, checkerboard harmonic ( $\xi = \pi$ ). The condition $G \le 1$ always holds.
In two and three dimensions the explicit scheme needs $r_x + r_y + r_z \le \tfrac12$ , the same condition as for a uniform grid in the lesson on the pressure equation.

For the explicit scheme $r \le \tfrac12$ means $\Delta t \le \Delta x^2/(2\eta)$ , the limit from the lesson Pressure equation for single-phase flow. Schemes with $\sigma \ge \tfrac12$ are stable at any step, but this does not mean that they damp fast harmonics well (more on this below, for the Crank–Nicolson scheme).

The analysis of harmonics links two views of stability from the earlier lessons. Each harmonic behaves like a separate “tank” with its own characteristic time $\tau_\xi = \Delta x^2/(4\eta\sin^2(\xi/2))$ . An explicit step multiplies it by $1 - \Delta t/\tau_\xi$ , exactly like the drawdown of one tank in the lesson on the simulation model, and it is stable for $\Delta t \le 2\tau_\xi$ . The “fastest” tank is the checkerboard harmonic with $\tau_\pi = \Delta x^2/(4\eta)$ . From this we get $\Delta t \le \Delta x^2/(2\eta)$ , the same limit that the checkerboard pattern gave in the lesson on the pressure equation.

CalculatedThe amplification factor over one step for harmonics of different wavelengths, for

r \approx 4

(cells of 250 m, a step of 30 days, hydraulic diffusivity 8412 m²/d), and the exact decay

\e^{-r\xi^2}

. The explicit scheme goes beyond the stability limit already for fairly long waves. The implicit scheme damps short waves, and damps long waves less than the exact solution does. In the Crank–Nicolson scheme short waves decay slowly and change sign.

Worked exampleAmplification factors of the three schemes on a grid with cells of 250 m

Given (assumed data)

Grid and rock (lesson Flow regimes): $\Delta x$ = 250 m; $\eta$ = 8412 m²/d
Time steps: 3, 6 and 30 d
Harmonics: a smooth one, $\xi = \pi/8$ (a wave of 16 cells = 4 km), and the checkerboard one, $\xi = \pi$

Solution

The numbers $r = \eta\Delta t/\Delta x^2$ : 0.4038, 0.8076 and 4.04. The limit of the explicit scheme is $r \le 0.5$ , that is $\Delta t \le \Delta x^2/(2\eta) = 3.71$ d.
The smooth harmonic, $\sin^2(\pi/16) = 0.0381$ . Step of 3 d: explicit 0.9385, implicit 0.9421, Crank–Nicolson 0.9404; exactly, $\e^{-r\xi^2} = 0.9396$ . All three are close to the exact value.
The checkerboard harmonic, $\sin^2(\pi/2) = 1$ . Step of 3 d: explicit $1 - 4\cdot0.4038 = -0.615$ , which decays while changing sign; implicit 0.382; Crank–Nicolson 0.106.
Step of 6 d: the smooth harmonic of the explicit scheme is 0.877 against the exact 0.883, still good. The checkerboard one is $1 - 4\cdot0.8076 = -2.23$ : in 10 steps it grows about 3000 times, and in 30 steps about 2.8·10¹⁰ times.
Step of 30 d, implicit scheme: the smooth one is 0.619 against the exact 0.536, and the checkerboard one is 0.058, almost fully damped.
Step of 30 d, Crank–Nicolson scheme: the smooth one is 0.530, almost exact. The checkerboard one is $(1 - 2\cdot4.04)/(1 + 2\cdot4.04) = -0.78$ : after 5 steps 0.29 of its initial amplitude remains (the sign changes at every step), while in the exact solution it disappears in the very first step.

Answer: On this grid the explicit scheme is stable for steps up to 3.71 d. At 6 d the checkerboard harmonic grows 2.23 times per step. At a step of 30 d the implicit scheme damps the checkerboard harmonic (0.058), and the Crank–Nicolson scheme hardly does (−0.78), but it reproduces the smooth harmonic more accurately (0.530 against 0.619, with the exact 0.536).

What it means. The factors explain the behaviour of the schemes without a calculation over cells. Instability is $|G| > 1$ for short waves. The lag of the implicit scheme is its $G$ for long waves being larger than the exact one: the disturbance decays more slowly than it should. The “ringing” of the Crank–Nicolson scheme is a $G$ close to −1 for short waves. The initial data and the round-off errors contain all harmonics, so an unstable harmonic will show up sooner or later, even if we do not see it in the first steps. The data in this example are assumed for illustration.

The analysis of harmonics does not take boundaries and heterogeneity into account. Usually they make the limit a little less strict: in the strip of the lesson Pressure equation for single-phase flow the smallest characteristic time of the cells is 3.3 d, and the exact limit is about 4.1 d. So in practice we use the local condition for the smallest and most conductive cells. It is reliable, and we check it for each cell separately.

The Crank–Nicolson scheme

For $\sigma = \tfrac12$ the flows during the step are taken as the average of the flows calculated from the pressures at the start and at the end of the step. In effect, we calculate the integral over time by the trapezoid rule (the average of the two ends). This is the Crank–Nicolson scheme, proposed in 1947 for the heat conduction equation. In the notation of the cell equation:

The Crank–Nicolson scheme for the cell equation

S_i\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \frac12\sum_j a_{ij}\left(p_j^{n+1} - p_i^{n+1}\right) + \frac12\sum_j a_{ij}\left(p_j^{n} - p_i^{n}\right) - q_i

where:

$S_i,\ a_{ij}$ storage capacity of the cell and the link coefficients with the neighbours, m³/MPa, m³/(d·MPa)
$q_i$ withdrawal by a well from the cell at reservoir conditions, m³/d
$p^{n},\ p^{n+1}$ pressures at the start and at the end of the step, MPa

Units: metric oilfield units, as in the lesson on the pressure equation

Holds when

A linear equation. The withdrawal $q_i$ is constant during the step.
The system for the pressures at the end of the step is tridiagonal (in one dimension), with the diagonal $S_i/\Delta t + \tfrac12\sum_j a_{ij}$ .

The work per step is the same as for the implicit scheme, one system of equations, but the error in time is of second order. The scheme is stable at any step: $|G| \le 1$ for all harmonics.

The weak point is visible in the amplification factor. For $\sigma = \tfrac12$ and the checkerboard harmonic, $G = (1 - 2r)/(1 + 2r)$ . With a large step ( $r \gg 1$ ) it tends to −1: short waves hardly decay and change sign at every step. While the solution is smooth, it has few short waves, and this goes unnoticed. But an abrupt change brings short waves into the solution. Examples are starting or shutting in a well, a change of rate, or a sharp initial pressure difference. Then decaying saw-tooth oscillations appear in the cells next to the place of the change. The implicit scheme damps such waves at once ( $G \to 0$ ). This property, the strong damping of the fastest components (those that change most quickly from cell to cell), is called L-stability.

Worked exampleOne step of 10 days with the Crank–Nicolson scheme on a strip of four cells

Given (assumed data)

Strip (lesson Pressure equation for single-phase flow): $S$ = 316.875 m³/MPa; $a_{12}$ = 68.24, $a_{23}$ = 27.296, $a_{34}$ = 28.433 m³/(d·MPa); withdrawal from cell 4 is 130 m³/d at reservoir conditions; initial pressure 20 MPa
Step: $\Delta t$ = 10 d, three times larger than the smallest characteristic time of the cells (3.3 d)

Solution

At the start of the step the pressure is uniform and there are no flows. So for the pressure changes $\delta_i$ the system is the same as for the implicit scheme, but with half the link coefficients: $\tfrac12 a$ = 34.120, 13.648 and 14.217. The diagonal is $b_i = S/\Delta t + \tfrac12\sum_j a_{ij}$ with $S/\Delta t = 31.6875$ : 65.8075, 79.4555, 59.5525 and 45.9045. The right-hand side is $(0;\ 0;\ 0;\ -130)$ .
Forward pass of the Thomas algorithm: $e_1 = 34.120/65.8075 = 0.5185$ ; $w_2 = 79.4555 - 34.120\cdot0.5185 = 61.76$ , $e_2 = 13.648/61.76 = 0.2210$ .
$w_3 = 59.5525 - 13.648\cdot0.2210 = 56.54$ , $e_3 = 14.217/56.54 = 0.2515$ ; $w_4 = 45.9045 - 14.217\cdot0.2515 = 42.33$ .
Backward pass: $\delta_4 = -130/42.33 = -3.071$ ; $\delta_3 = 0.2515\cdot\delta_4 = -0.772$ ; $\delta_2 = 0.2210\cdot\delta_3 = -0.171$ ; $\delta_1 = 0.5185\cdot\delta_2 = -0.089$ MPa.
Pressures after 10 days: 19.91, 19.83, 19.23 and 16.93 MPa. Exact solution of the cell equations: 19.92, 19.82, 19.11 and 17.04 MPa. The implicit scheme with the same step gave 17.39 MPa in cell 4.
Step of 20 d: the pressure in cell 3 from the Crank–Nicolson scheme deviates from the exact one alternately in opposite directions, by +0.24, −0.15 and +0.03 MPa after 20, 40 and 60 days. For the implicit scheme the deviations have one sign: +0.06, +0.11 and +0.08 MPa.

Answer: After 10 days, in the cell with the well, the Crank–Nicolson scheme gives 16.93 MPa against the exact 17.04 MPa: an error of −0.11 MPa, against +0.35 MPa for the implicit scheme with the same step.

What it means. The second order in time shows in the cell with the well: for the same work, one tridiagonal system per step, the error there is three times smaller. In the neighbouring cell 3 the picture is the opposite. Already after 10 days the Crank–Nicolson scheme is wrong by more than the implicit scheme (+0.12 against −0.06 MPa). With a step of 20 days the pressure in cell 3 “swings” around the exact solution, and its first deviations are larger than for the implicit scheme. Starting a well is an abrupt change, and the fast harmonics of this scheme decay slowly and change sign. The implicit scheme has a larger error in the cell with the well, but its errors have one sign and no oscillations. For forecasts in which the operating modes of the wells change often, this matters more. The data in this example are assumed for illustration.

CalculatedA test problem: in a closed strip 2 km long, the pressure in the left half is 2 MPa higher than in the right half. The profile after 10 days, calculated in one step of 10 days on a grid of 20 cells of 100 m (

r \approx 8.4

). The implicit scheme gives a smooth profile. The Crank–Nicolson scheme “flips” the pressure difference in the two middle cells, and at the next step it will flip back.

There are several ways to remove the oscillations. The most reliable is to make two or three steps with the implicit scheme after every abrupt change, and only then to go over to $\sigma = \tfrac12$ : the implicit steps damp the short waves. Reducing the step after a change in how the wells are run also helps. A weight a little above one half (for example, 0.55) weakens the “ringing”: with a large step the checkerboard harmonic is multiplied not by −1 but by about $-(1 - \sigma)/\sigma \approx -0.82$ . It also keeps the error in time small, but it does not remove the oscillations completely. Reservoir simulators usually use the first-order implicit scheme in time. The multiphase equations are nonlinear, the operating modes of the wells change often, and reliability matters more than second order. They reach the needed accuracy by the choice of the step (IMPES and fully implicit schemes).

A convergence study for the grid and for the step

The lesson on grids describes the check of convergence in the grid: calculations on several grids, each twice as fine as the one before, the observed order (the order that the results actually show), and extrapolation to an infinitely fine grid. The theory of finite-difference schemes adds two things to it: what to expect from the observed order, and how not to confuse the grid error with the time-step error. For a smooth solution the error of a result is the sum of two parts:

Error model for a convergence study

e = R - R_{\text{exact}} \approx C_x\,\Delta x^{p} + C_t\,\Delta t^{q}

where:

$R,\ R_{\text{exact}}$ result of the calculation (pressure, rate, breakthrough time) and its exact value, in units of the result
$p,\ q$ orders in space and in time: for the pressure equation $p$ = 2; $q$ = 1 for the implicit scheme and 2 for the Crank–Nicolson scheme, —
$C_x,\ C_t$ coefficients that depend on the solution (its derivatives) and on the scheme, but not on the steps, —

Units: in units of the result

Holds when

The steps are small enough for the leading terms of the error to dominate (the asymptotic range).
The solution is smooth. Near wells, at displacement fronts and at sharp boundaries of properties, the observed order is lower than the theoretical one.
The errors of the solution of the equations at each step (the tolerances of the iterations) are much smaller than the error of the scheme.

If we refine only the grid, the error tends not to zero but to $C_t\Delta t^q$ . If we refine only the step, it tends to $C_x\Delta x^p$ . So we check the step and the grid separately, or refine them together in a consistent way.

The grid at a fixed step. The differences of the results on three grids still show the order $p$ , because the constant part $C_t\Delta t^q$ cancels in them. But the extrapolation gives the answer for this step, not the exact one.
The step on a fixed grid. The error falls to a “floor”, the grid error. It is useless to reduce the step further.
Consistent refinement. If the grid and the step are made finer together, with $\Delta t \propto \Delta x$ , the implicit scheme finally converges with first order, because the error in time falls more slowly. To keep second order, we need a step $\Delta t \propto \Delta x^2$ (as for the explicit scheme, but now for accuracy and not for stability) or a scheme of second order in time.

How to check the convergence of a calculation

1
Choose the indicators
Take the quantities for which the model is built: the bottomhole pressure, the cumulative production by a date, the water breakthrough time. Different indicators converge differently: pressures usually converge faster than the breakthrough time.
2
Check the time step
On the working grid, run calculations with steps $\Delta t$ , $\Delta t/2$ and $\Delta t/4$ . The differences of the results show the order and the size of the step error. For first order, the error of the calculation with step $\Delta t/2$ is about equal to the difference between its result and the result with step $\Delta t$ .
3
Check the grid
Use three grids, each twice as fine as the one before, either with a step whose error is already small or with a step reduced together with the cells. The observed order from the three grids and the extrapolation follow the formula of the lesson Grids in reservoir simulation models.
4
Compare the order with the theory
An observed order that differs noticeably from the theoretical one means that the grids are still coarse (not in the asymptotic range), the solution is not smooth, or there is a mistake in the calculation, for example in the transmissibilities or the wells.
5
Compare the error with the purpose
We compare the estimated error with the accuracy that the decision needs: with the difference between development options, and with the uncertainty of the input data.

CalculatedThe error of the amplitude in the test problem (after 30 days) against the time step. Without the grid error, the implicit scheme converges with first order (slope 1), and the Crank–Nicolson scheme with second order (slope 2). On a grid of 8 cells, reducing the step runs into the grid error, 4.3 kPa.

Terms and symbols

Finite-difference method: Replacing derivatives by differences of values at grid nodes. It turns a differential equation into a system of algebraic ones, a finite-difference scheme.
Stencil: The set of grid nodes (in space and time) that enter the difference formula of one equation.
Truncation error $\psi$: The residual that the exact solution gives when it is substituted into the difference equation. Its leading term is proportional to a power of the step.
Order of accuracy: The exponent of the step in the leading term of the error. When the step is halved, the error falls about $2^m$ times.
Stability: The scheme does not amplify errors from step to step. Together with approximation it ensures convergence (the Lax equivalence theorem).
Amplification factor $G$: How many times the scheme changes the amplitude of a harmonic over a step. Stability means $|G| \le 1$ for all harmonics.
Mesh Fourier number $r$: $\eta\Delta t/\Delta x^2$ . The explicit scheme of the pressure equation is stable for $r \le 1/2$ .
Weighted scheme: The right-hand side is a weighted average of the values at the start and at the end of the step, with the weight $\sigma$ at the end: 0 is explicit, 1 is implicit, ½ is Crank–Nicolson.
Crank–Nicolson scheme: A scheme with weight ½: second order in time, stable at any step, but with a large step it damps fast components poorly.

Common mistakes and misconceptions

Wrong

“A scheme of second order is always more accurate than a scheme of first order.”

In fact

The order tells how fast the error falls as the step gets smaller. With a large step the Crank–Nicolson scheme gives oscillations after a well is started, which the implicit scheme does not have, and at displacement fronts the order of any scheme falls.

Wrong

“The first tens of steps look reasonable, so the scheme is stable.”

In fact

An unstable harmonic grows out of round-off errors: with a factor of −2.2 it grows three thousand times in 10 steps and tens of billions of times in 30. Meanwhile the smooth part of the solution may be calculated correctly right up to the “explosion”. We check the step against the stability condition.

Wrong

“The result stopped changing when we refined the grid, so it is exact.”

In fact

At a fixed step, the grids converge to the solution with that step. In the test problem of the lesson, the extrapolation over three grids gave second order and was still wrong by 0.010 MPa, which is the error of the 3-day step. We check the step separately.

Wrong

“The iterations converged at all steps, so the model has converged.”

In fact

The convergence of the iterations means that the difference equations are solved. How far these equations are from the differential ones is shown only by refining the grid and the step.

Wrong

“The smaller the step of numerical differentiation, the more accurate the derivative.”

In fact

Only for exact values. When the values have an error $\varepsilon$ , an error of the order of $\varepsilon/\Delta x$ is added to the error of the formula, and with a very small step the noise dominates. That is why the derivative of measured pressure is calculated from points spaced apart.

In depth: advanced detail

Derivation of the truncation error of the weighted scheme

We substitute the exact solution into the scheme and expand all the values in a Taylor series around the point $(x_i,\ t^{n+1/2})$ , the middle of the step. The difference in time is central with respect to the middle of the step, so $(p^{n+1} - p^{n})/\Delta t = \partial_t p + \tfrac{\Delta t^2}{24}\,\partial_t^3 p + \ldots$ At each moment the second difference equals $\partial_x^2 p + \tfrac{\Delta x^2}{12}\,\partial_x^4 p + \ldots$ , and the values at the start and at the end of the step differ from the values in the middle by $\mp\tfrac{\Delta t}{2}\,\partial_t$ plus terms of the order of $\Delta t^2$ .

Expansion of the weighted right-hand side

\sigma\,\delta^2 p^{n+1} + (1 - \sigma)\,\delta^2 p^{n} = \partial_x^2 p + \left(\sigma - \tfrac12\right)\Delta t\,\partial_t\partial_x^2 p + \frac{\Delta x^2}{12}\,\partial_x^4 p + O\!\left(\Delta t^2 + \Delta x^2\Delta t + \Delta x^4\right)

where:

$\delta^2 p$ three-point second difference divided by $\Delta x^2$ , MPa/m²
$\partial_t,\ \partial_x$ partial derivatives with respect to time and the coordinate in the middle of the step, —

Units: any consistent units

Holds when

A smooth solution, constant coefficients, a uniform grid.

We multiply by $\eta$ and subtract from the difference in time. The exact solution satisfies $\partial_t p = \eta\,\partial_x^2 p$ , so the leading terms cancel, and $\eta\,\partial_t\partial_x^2 p = \partial_t^2 p$ . What remains is $\psi = (\tfrac12 - \sigma)\Delta t\,\partial_t^2 p - \tfrac{\eta\Delta x^2}{12}\partial_x^4 p$ , the formula of the main section.

A non-uniform grid

Let the distances from node $i$ to its right and left neighbours be $h_+$ and $h_-$ . Taylor expansions again give the generalisation of the three-point formula and its error:

The second difference on a non-uniform grid of nodes

\frac{2}{h_+ + h_-}\left(\frac{p_{i+1} - p_i}{h_+} - \frac{p_i - p_{i-1}}{h_-}\right) = p''_i + \frac{h_+ - h_-}{3}\,p'''_i + O(h^2)

where:

$h_+,\ h_-$ distances to the right and left neighbouring nodes, m
$h$ characteristic step of the grid, m

Units: any consistent units

Holds when

A grid of nodes (values at points). For a block-centred grid the analysis is more complicated, but the practical conclusion is the same.
A smooth function.

If neighbouring steps differ by a fixed factor, the leading term is proportional to $h$ and the order falls to first. If the step changes smoothly, $h_+ - h_-$ is of the order of $h^2$ , and second order is kept. Hence the rule: we change the sizes of neighbouring cells gradually and not abruptly (local grid refinement).

A scheme that converges to the wrong equation

We cannot take approximation for granted. A classic example is the DuFort–Frankel scheme. Take an explicit scheme with a central difference in time over two steps, $(p_i^{n+1} - p_i^{n-1})/(2\Delta t)$ . For the pressure equation it is unstable at any step. If in the second difference we replace $2p_i^n$ by $p_i^{n+1} + p_i^{n-1}$ , the scheme stays explicit and becomes stable at any step. But the Taylor expansion adds to its error a term $\eta\,(\Delta t/\Delta x)^2\,\partial_t^2 p$ . It vanishes only if the time step is reduced faster than the grid step. If we refine them proportionally, $\Delta t/\Delta x = \beta$ , the scheme converges to the solution of a different equation, $\partial_t p + \eta\beta^2\partial_t^2 p = \eta\,\partial_x^2 p$ . This equation has an “inertial” term that does not exist in flow through porous rock. The conclusion: we check approximation for the way of refinement that is really used.

A two-dimensional grid: the five-point stencil and the numbering of cells

In two dimensions the second differences in $x$ and in $y$ add up, and the equation of cell $(i, j)$ contains the cell itself and its four neighbours that share a face: the five-point stencil. For the explicit scheme the analysis of harmonics gives $G = 1 - 4r_x\sin^2(\xi/2) - 4r_y\sin^2(\zeta/2)$ , with phases $\xi$ and $\zeta$ along the two axes. The worst case is the checkerboard harmonic along both axes, from which $r_x + r_y \le \tfrac12$ .

To solve the implicit system, the cells are numbered with a single index. With the natural row-by-row numbering, cell $(i, j)$ gets the number $m = (j - 1)N_x + i$ , where $N_x$ is the number of cells along $x$ . The neighbours in $x$ have the numbers $m \pm 1$ , and the neighbours in $y$ have $m \pm N_x$ . The matrix of the system has five non-zero diagonals: the main one, the two next to it, and two that are $N_x$ places away. The width of this “band” is set by the number of cells in a row. So for direct elimination methods (they solve the system by elimination in a fixed number of operations) it pays to number the cells first along the short side of the model. For a model of 150 × 30 cells the band narrows from 150 to 30, and the amount of work of banded elimination, which is proportional to the square of the band width, falls about 25 times. How such systems with millions of unknowns are solved is in the lesson Solving systems of linear equations in simulators.

Connections to other topics

Pressure equation for single-phase flow — the cell equations, the explicit and implicit schemes, the largest step and the Thomas algorithm.
What a reservoir simulation model is and why it is needed — the explicit and implicit step for one “tank”, the simplest case of the analysis of harmonics.
Grids in reservoir simulation models — the check of convergence in the grid and the extrapolation to an infinitely fine grid; numerical diffusion of the upstream scheme.
Diffusivity equation — the diffusivity equation as the model equation for finite-difference schemes.
Two-phase flow equations — the pressure and saturation equations to which the same schemes are applied.
IMPES and fully implicit schemes — schemes for multiphase nonlinear equations, Newton iterations and the choice of the step.
Solving systems of linear equations in simulators — how the sparse systems that implicit schemes give are solved.
Numerical dispersion — the error of schemes for the transport of saturation, and the smearing of fronts.
Diagnosis with the pressure derivative — numerical differentiation of noisy measurements.

Check yourself

1
A central difference with a step of 40 m is wrong by 0.6%. What error should we expect for a step of 10 m (a smooth function, exact values)?
1. A.0.6%
2. B.0.15%
3. C.0.04%
4. D.0.01%
Show the answer
Correct answer: C. 0.04%
The central difference is of second order: the error is proportional to the square of the step. The step fell 4 times, so the error fell 16 times: $0.6/16 \approx 0.04$ %. The option 0.15% would correspond to first order.
2
A scheme with weight $\sigma$ = 0.25 solves the pressure equation on a grid with cells of 100 m, with a hydraulic diffusivity of 20,000 m²/d. At what largest step is it stable, and how many times larger is this step than that of the explicit scheme?
Show the answer
The condition for $0 < \sigma < \tfrac12$ is $r \le 1/(2(1 - 2\sigma)) = 1/(2\cdot0.5) = 1$ , that is $\Delta t \le \Delta x^2/\eta = 10{,}000/20{,}000 = 0.5$ d. For the explicit scheme $r \le 0.5$ , so $\Delta t \le 0.25$ d. The limit is twice as large, but the work per step is larger too: we have to solve a system of equations, as in the implicit scheme.
3
Which scheme, with a large step after a well is started, gives decaying pressure oscillations in the cells next to the well?
1. A.Implicit
2. B.Crank–Nicolson
3. C.Explicit with $r < 1/2$
4. D.None: all stable schemes are monotonic
Show the answer
Correct answer: B. Crank–Nicolson
For the Crank–Nicolson scheme the amplification factor of the shortest harmonic, $(1 - 2r)/(1 + 2r)$ , is close to −1 for $r \gg 1$ : it hardly decays and changes sign at every step. Starting a well creates an abrupt change, rich in short harmonics. For the implicit scheme the factor $1/(1 + 4r)$ is positive and small. The explicit scheme with $r < 1/2$ is not used with a large step. Its new cell pressure is a weighted average of the pressures of the cell itself and its neighbours, with non-negative weights (Pressure equation for single-phase flow), so the pressure changes monotonically.
4
The bottomhole pressure was calculated on three grids, each twice as fine as the one before, with the same step of 5 d: 14.62, 14.95 and 15.03 MPa. What is the observed order, what is the extrapolated value, and what else must be checked?
Show the answer
The changes are 0.33 and 0.08 MPa, and the order is $\log_2(0.33/0.08) = \log_2 4.1 \approx 2.0$ , as for the second difference. The extrapolation is $15.03 + 0.08/3 \approx 15.06$ MPa. But this is the value for a step of 5 d: the grids do not see the step error. We need calculations with a step of 2.5 d (and 1.25 d) on one of the grids to estimate it.
5
A scheme approximates the linear pressure equation with the error $O(\Delta t + \Delta x^2)$ and is stable. What can we say about its solution as the grid and the step get smaller?
1. A.It converges to the exact one with first order in time and second order in space
2. B.It converges only for $r \le 1/2$
3. C.It converges, but the order cannot be predicted
4. D.It may not converge: another check is needed
Show the answer
Correct answer: A. It converges to the exact one with first order in time and second order in space
By the Lax equivalence theorem, for a linear problem approximation and stability give convergence. The error of the solution is bounded by the calculation time multiplied by the largest truncation error (up to a constant factor), so it falls with the same orders. The condition $r \le 1/2$ belongs only to the explicit scheme. Here the question already says that the scheme is stable.
6
Why may a convergence study only in the grid, with a fixed time step, be insufficient even if the observed order is “right”?
Show the answer
The error is the sum of a spatial part $C_x\Delta x^p$ and a temporal part $C_t\Delta t^q$ . At a fixed step the temporal part is constant. It cancels in the differences of the results on different grids, so the order comes out right, and the extrapolation converges to the solution with this step. In the test problem of the lesson an error of 0.010 MPa remained in this way, almost 2% of the amplitude. We check the step with a separate series of calculations.

Updated: 10 October 2026