A different equation on every grid

The stencil that is not symmetric established the convection field’s organising result: upwinding is not a better approximation of the problem posed, it is an exact discretisation of a different one. The upwind matrix at diffusion ε is the central-difference matrix at diffusion ε(1 + Pe), entry for entry, and the upwind answer sits on the exact solution of that thicker-layered problem rather than near the exact solution of the thin-layered one it was asked about.

What that leaves unsaid is which different problem. The phrase invites a reader to imagine one substitute equation — the equation upwinding really solves — sitting beside the equation that was posed, at a fixed distance from it, waiting to be named. There is no such equation. What upwinding is exact for is a sequence of equations, one for every mesh, and the index of that sequence is h.

The added diffusion is εPe, and Pe is h/2ε, so the ε cancels and what is left is h/2. That is one line of algebra and it is stated in the essay above. Stating it and showing it are different things, and with a single grid it cannot be shown: at n = 31 the quantity h/2 is the number 0.015625, and a number that never moves is indistinguishable from a constant somebody chose. It becomes a variable only when the mesh is swept, and only then does the claim have any content — that ε, the physical quantity the whole problem is about, is absent from the perturbation the scheme commits.

The diffusion asked for is not in the diffusion added #

Hold the mesh at 31 points and move ε across two and a half decades. The cell Péclet number goes 0.07813, 0.1563, 0.3125, 0.7813, 3.125, 7.813, 15.63 at ε = 0.2, 0.1, 0.05, 0.02, 0.005, 0.002 and 0.001 — a factor of exactly two hundred from end to end, and the interval crosses the threshold at which the central-difference matrix stops being an M-matrix.

The artificial diffusion εPe across those same seven problems reads 0.015625 at every one of them. Not 0.0156 and 0.0157; the same double, to the last bit, seven times. It is h/2 on that grid, and h/2 knows nothing about ε.

That is worth separating from the thing it is often confused with, because a reader who has just met the identity will reasonably expect the perturbation to grow as the problem gets harder, and it does not. What grows is the ratio. The factor (1 + Pe) by which the effective diffusion exceeds the one asked for reads 1.078, 1.156, 1.313, 1.781, 4.125, 8.813 and 16.625 across the same seven values. At ε = 0.2 the scheme is solving a problem 7.8% more diffusive than the one posed; at ε = 0.001 it is solving one sixteen and a half times more diffusive. The absolute perturbation is identical in all seven cases and the relative one spans a factor of fifteen.

So ε does not decide how large the perturbation is. It decides how large the perturbation looks against the physics, which is a different quantity and is the one a reader cares about — and it is also the reason the identity feels like a statement about hard problems when it is a statement about coarse grids.

A sequence indexed by h, and nothing else #

Now fix ε at 0.005 and move the mesh instead. Six grids, from 15 interior points to 511, each a doubling of the last:

Two columns of that table are the same column. The added diffusion and h/2 agree not to a tolerance but as floating-point numbers: their ratio minus one is exactly zero at all six meshes, which is the strongest form the claim can take and is stronger than the norm comparison beside it.

The last column is the identity itself, measured at each mesh rather than argued. The relative Frobenius difference between the upwind operator and the central operator for the perturbed equation runs 2.57·10⁻¹⁷, 1.26·10⁻¹⁶, 2.59·10⁻¹⁷, 3.84·10⁻¹⁷, 0 and 0. There is no trend in that column and there should not be, because the residue is the rounding committed in forming ε(1 + Pe) and multiplying it by h⁻², not an error in the identity. Two of the six meshes give matrices that agree bit for bit. The identity is not asymptotic, not approximate and not better on fine grids: it holds wherever it is checked.

The identity is also stronger than an identity between two matrices, which is worth checking because a linear system has a right-hand side and the two schemes do not share one. The boundary condition u(1) = 1 enters the last equation differently in each: the central scheme contributes ε/h² − 1/2h and the upwind scheme contributes ε/h². Evaluate the first at the perturbed diffusion and the −1/2h cancels against the +1/2h that the extra h/2 of diffusion brings, leaving exactly ε/h². So the two right-hand sides agree as well, at a relative difference of 5.20·10⁻¹⁶, 1.73·10⁻¹⁶, 1.73·10⁻¹⁶, 1.73·10⁻¹⁶, 0 and 0 across the six meshes — one unit in the last place at four of them and nothing at the other two. The upwind scheme and the central scheme at ε + h/2 are the same linear system, not merely the same matrix, and their computed solutions differ by at most 7.5·10⁻¹⁶ and by exactly zero on the two finest grids.

That closes off the obvious way the claim could have been true and useless. A matrix identity with different right-hand sides would give two different answers to two different problems and would say nothing about which equation the upwind answer belongs to; the systems being equal is what makes the substitution total.

And a further reading of the same table, which is the one that makes “a sequence” the right word. No member of the sequence is the equation that was posed. The added diffusion is h/2 and h is positive at every finite mesh, so ε + h/2 is strictly larger than ε on every grid a code can build. The sequence converges to the posed equation and never contains it, which is a different situation from a scheme that is exact for the right problem at some resolution and wrong at others. There is no grid at which upwinding solves the problem asked about.

The distance to the posed problem on that grid is worth reading twice. A relative error of nineteen thousand is not a bad approximation; it is not an approximation at all. The exact solution at ε = 0.005 is essentially zero across almost the whole interval and rises inside a layer of width 0.005, and the upwind answer rises inside a layer seven times wider, so the ratio between them where the first is small is enormous. That number is what the identity predicts and no solver produces.

The mesh crosses the threshold from the other side #

There is a point in that sequence where something visible happens to the other scheme. The cell Péclet number falls with h, so refining walks the problem back across the threshold at Pe = 1 — the point at which the upper off-diagonal −ε/h² + 1/2h changes sign, the matrix stops being an M-matrix and the central-difference answer stops obeying the maximum principle the equation obeys. At ε = 0.005 that crossing sits between 63 points, where Pe is 1.563, and 127, where it is 0.7813.

The crossing matters here because it is where the two repairs stop being alternatives. Above the threshold there is a decision to make: accept an answer that leaves [0, 1] at half the grid points, or accept a different equation. Below it there is none, because central differencing is monotone and second-order accurate and gives up nothing. The whole reason upwinding exists is the region above the line, and refining the mesh removes the region rather than improving the answer inside it.

That figure is drawn at the mesh the sequence has reached rather than at the ε the earlier essay swept, and the two readings are different claims. Sweeping ε at a fixed grid shows the threshold as a property of the problem; sweeping the grid at a fixed ε shows it as a property of the discretisation, and only the second is what the Péclet number actually is. The refusal this essay publishes is exactly that: an assertion that Pe does not depend on the mesh is fed 6.25 on 15 points against 0.390625 on 255, at the same ε, and is required to fail.

The equation upwinding solves is always on the safe side of it #

There is a consequence of the identity that has to be true and is not usually said, and it explains the one property upwinding is chosen for.

The upwind scheme is monotone at every mesh and every ε. The standard argument for that is a statement about the matrix: both off-diagonals are negative whatever h and ε are — at 31 points and ε = 0.005 they are −37.12 and −5.12, against −21.12 and +10.88 for central differences on the same problem — so the matrix is an M-matrix, its inverse is entrywise non-negative, and a discrete maximum principle follows.

The identity gives the same fact a second derivation with nothing about matrices in it. The equation upwinding is exact for has diffusion ε + h/2, so its own cell Péclet number is

Pe′  =  h / 2(ε + h/2)  =  (h/2) / (ε + h/2)

which is a positive number divided by itself plus ε, and is therefore strictly less than one for every ε > 0 and every h. Measured along the mesh sweep at ε = 0.005 it reads 0.862, 0.758, 0.610, 0.439, 0.281 and 0.163 at n = 15 to 511; measured along the ε sweep at 31 points it reads 0.0725, 0.135, 0.238, 0.439, 0.758, 0.887 and 0.940 at ε = 0.2 down to 0.001. It approaches one from below as ε falls and never reaches it: at ε = 10⁻¹² on that grid it is 0.99999999994.

So upwinding is central differencing applied to a problem that has been moved just far enough to be inside the region where central differencing is safe. That is two routes to a number applied to the maximum principle — one through the sign of an off-diagonal entry, one through a threshold on the equation — and the second says something the first does not. The threshold is never crossed with room to spare. At ε = 0.001 on 31 points the modified problem sits at Pe′ = 0.940, which is under the line by six per cent, and the scheme’s guarantee is that it is under the line rather than that it is comfortably under it.

What converges at first order is the equation #

The convection field has already fitted the two schemes’ orders of accuracy, at an ε large enough that the boundary layer is resolved on every grid: central differences at 2.011 and upwinding at 0.988, over the five meshes from 15 points to 255. Second order and first order, fitted rather than quoted, and the standard reading is that upwinding is the less accurate discretisation.

That reading is wrong, and the identity is what makes it wrong. If the upwind operator is the central operator for the equation with diffusion ε + h/2, then it is a second-order discretisation — of that equation. Its first-order behaviour has to be coming from somewhere else. The measurement that locates it splits the upwind error into two pieces which can be computed separately.

The first piece is the discretisation’s own error: the distance from the upwind answer to the exact solution of the equation the upwind operator is exact for. The second is a modelling error with no discretisation in it at all: the distance between the exact solution at ε + h/2 and the exact solution at ε, both in closed form, sampled on the same grid points. At ε = 0.5 the four fitted orders read

central differences, against the posed problem        2.011
upwinding, against the equation it discretises        1.958
upwinding, against the posed problem                  0.988
the two exact solutions, against each other           0.994

Upwinding is second order. Against its own equation it converges at 1.958, which is the same rate central differencing achieves against the posed one, and the two schemes are the same scheme. What converges at first order is the last line — a comparison between two exact solutions of two continuous problems, in which no matrix appears.

The shares settle it. Across the five meshes the modelling error alone accounts for 1.0190, 1.0099, 1.0051, 1.0026 and 1.0013 times the total upwind error, so by 255 points the gap between the two equations explains the whole of the error to within 0.13%, and the discretisation’s own contribution has fallen out of sight. At ε = 0.1 the same four fits read 2.021, 1.881, 0.960 and 0.985, and the shares run 1.0814 down to 1.0064: the same decomposition, arriving more slowly because the layer is thinner.

So the familiar sentence — upwinding is first-order accurate — is not a statement about how faithfully the scheme represents its equation. It is a statement about how fast the family of equations being solved approaches the one that was asked for, and the rate is first order because the family is indexed by h/2.

This is the exact answer to a nearby problem read one level further out than it is usually read, and read with the perturbation named rather than bounded. Backward error is normally a statement about an algorithm: a computed answer is described by the problem it is the exact answer to, and the size of the change is what the algorithm is answerable for. Here the object with the backward error is the discretisation, the perturbed problem is a differential equation rather than a matrix, and the perturbation is not merely small but written down: ε becomes ε + h/2. That is the same move an accuracy that is a backward error makes for a compression tolerance, where the knob a caller sets is the backward error itself rather than a bound on it.

It also explains the asymmetry between the two repairs that the earlier essay recorded and did not account for. Refining is the repair that always works because h is in the equation. Every other lever a numerical code has — a better solver, a tighter tolerance, a longer mantissa — acts on the arithmetic, and the arithmetic is not where the error is: the oscillating central answer solves its own system to 4.6·10⁻¹⁸ and the upwind answer solves its own to the same order. A small residual is not a small error separates an algorithm’s contribution from a problem’s; this separates both of them from a third thing upstream, which no quantity computed from the assembled matrix can see. Refining is the only lever that reaches it, and it reaches it because the discretisation parameter is a coefficient of the equation being solved.

Where the reading stops #

Three limits, each of which is a different essay’s subject.

The perturbation is h/2 for this upwinding and not for artificial diffusion in general. The diffusion that makes the answer exact adds ε·ξ·Pe with ξ = coth(Pe) − 1/Pe, which is a fraction of the same h/2 and depends on Pe and therefore on ε again. The clean statement here — a perturbation that contains no physics — is a property of the crudest scheme in the family, and buying accuracy back costs it.

It is a statement about one space dimension and one flow direction. In two dimensions the added diffusion is a tensor rather than a number, and the direction the diffusion does not go measures what a rank-one artificial diffusion leaves undamped across the flow. Exact along one axis measures what happens when the flow is not aligned with the grid, where the nodal exactness of the tuned scheme survives exactly one angle. Neither result is contained in the one-dimensional identity.

And refining is affordable here and often is not. At ε = 0.005 the mesh that brings the added diffusion below a tenth of ε is 999 points, which is nothing. At ε = 10⁻⁶, an ordinary physical value, Pe = 1 needs 499,999 points in one dimension and a tenth of ε needs 4,999,999 — before any second dimension is added. The repair that always works is the repair that is frequently out of reach, which is the entire reason a scheme with a backward error exists.

What follows for a code that solves one of these #

A convergence study measures the equations, not the scheme. A table of errors against h on an upwind discretisation is reporting the rate at which ε + h/2 approaches ε. It is a real measurement and it is not a measurement of discretisation fidelity, and the two are usually reported under the same heading.

A tolerance is the wrong instrument for this error. Nothing in the solve is inaccurate, so tightening an iterative tolerance, switching to a direct method or raising the working precision moves none of it. The quantity to report beside the answer is h/2 against ε — one division — and it says how different the equation solved was from the equation posed.

The perturbed problem is inside the class asked about, which is not guaranteed. ε + h/2 is a convection–diffusion problem of exactly the posed kind, with a larger coefficient, so every intuition about the original transfers. A nearby problem of the wrong kind is the case where that fails — a backward error whose explaining perturbation leaves the structure the problem was defined by — and it is worth knowing that the convection case is the benign one rather than assuming it is the general one.

The model’s own error can dominate the model’s condition number. The condition number of the model measures how much a κ computed over the wrong class of perturbations over-states the difficulty of a structured problem. Here the analogue is sharper: the perturbation is not merely of a restricted class, it is known exactly, so the error it causes can be computed rather than bounded — which is what the fourth fitted order above is.

And a rate derived on the Laplacian is a rate for a different operator. A rate that is known in advance checks Jacobi, Gauss–Seidel and optimal SOR against closed forms belonging to a symmetric operator, and every one of those forms is derived for the equation posed rather than the equation upwinding solves. On this operator the diffusion coefficient in the closed form is ε + h/2, so a rate computed from ε is a rate for a matrix nobody assembled.

The general shape is one this collection has met in another field entirely. An error committed before the arithmetic is about a choice made when a problem is set up, which fixes the accuracy of everything downstream and cannot be revised afterwards. Choosing a mesh is that choice for a differential equation, and the reason it is easy to miss is that it does not look like a modelling decision. It looks like a resolution setting. It is a coefficient.