The diffusion that makes the answer exact

The convection essay leaves the field with two schemes and a threshold between them. Central differencing solves its own system to 4.6·10⁻¹⁸ and, past a cell Péclet number of exactly one, produces a solution that leaves the interval [0, 1] at sixteen of thirty-one points — an answer the differential equation forbids, computed correctly. Upwinding never leaves the interval, and it is solving a different problem: ε(1 + Pe) instead of ε, entry for entry to 1.3·10⁻¹⁶, with the added diffusion exactly h/2 and independent of ε.

Both of those are choices about a single number, and neither choice was made by asking what the answer was. Upwinding adds h/2 because that is what a one-sided difference happens to add. Central differencing adds nothing because that is what a symmetric difference happens to add.

This essay is about the scheme that picks the number on purpose, and about what its success is a property of.

The tuning is an exact statement, so it is checked at every stop the slider has rather than at the one the figure opens on.

Two things to watch as ε falls: the tuned column, which should not move, and the order of the other two, which does.

A third of a per cent of extra diffusion at the first stop and ten per cent at this one, and in both cases the nodal values come out exact. The next stop is the last one at which central differences are still the sensible ordinary choice:

Eleven nodes outside the range the solution is known to live in, on a scheme that was the more accurate of the two one stop earlier. Two more decades of convection make that worse for one scheme and better for the other:

The tuned scheme is exact at every stop. 2.78·10⁻¹⁶, 5.55·10⁻¹⁷, 5.55·10⁻¹⁷, 2.78·10⁻¹⁷, 6.94·10⁻¹⁸, 2.36·10⁻¹⁷, 5.80·10⁻¹⁷ — seven readings, all at the unit roundoff, across a hundredfold range of ε and a hundredfold range of Péclet number. Nothing about the exactness is asymptotic and nothing about it degrades; the nodal values are right to the last bit whatever the flow does.

ξ/Pe is exactly a third at small Péclet and falls away. 0.333, 0.333, 0.331, 0.321, 0.289, 0.219, 0.112 — the small-argument expansion of coth(Pe) − 1/Pe is Pe/3, and the first three stops confirm it to three digits. So the tuning is not a large correction where it is not needed: at Pe = 0.078 it adds 2.6% of the diffusion already there.

And two thresholds coincide at Pe = 1. Central differences are the more accurate of the two ordinary schemes up to Pe = 0.781 — 0.0868 against 0.1806 — and the worse of them from Pe = 1.563 — 0.2634 against 0.1985. The nodes outside [0, 1] appear at exactly the same step: 0 at 0.781 and eleven at 1.563. The accuracy crossover and the monotonicity threshold are the same threshold, which is not obvious and is what makes Pe = 1 a real boundary rather than a convention.

Upwind is non-monotone in the flow, which is the last surprise here. Its error runs 0.0254, 0.0506, 0.0922, 0.1806, 0.1985, 0.1360, 0.0602 — rising to a peak at Pe = 1.563 and then improving by a factor of three as the problem becomes more convection-dominated. The scheme is worst where the two effects are comparable and better at both extremes, so a code benchmarked at one Péclet number has learnt very little about it.

One formula #

ξ(Pe)  =  coth(Pe) − 1/Pe        ε̃  =  ε (1 + ξ·Pe)

Discretise with central differences at ε̃ rather than ε, and the computed solution equals the exact solution at every grid point.

Not “second order”, not “uniformly convergent in ε”, not “accurate for practical purposes”: exact, at the nodes, at every Péclet number tried, with the error at the level of the arithmetic that computed it.

The site has one other family with this property — the Hilbert matrix, whose inverse exact.js computes in BigInt rationals, and the Kac–Murdock–Szegő matrix, whose inverse is tridiagonal in closed form. Both are exact ground truths because somebody wrote down a formula for an answer. This one is a formula for a scheme, and the exactness is of the discretisation rather than of the solve.

Why it works, which is also why it is fragile #

The exact solution of −εu″ + u′ = 0 with u(0) = 0, u(1) = 1 is

u(x)  =  (e^(x/ε) − 1) / (e^(1/ε) − 1)

which, sampled at equally spaced points, is a geometric sequence plus a constant: uᵢ = (rⁱ − 1)/(rⁿ⁺¹ − 1) with r = e^(h/ε).

A three-point difference operator applied to a geometric sequence gives that sequence back multiplied by a quadratic in its ratio. So the discrete equation is satisfied exactly if the operator’s quadratic has a root at exactly r — and ξ is the value of the diffusion coefficient that puts it there. The derivation is two lines and it consists entirely of matching the discrete ratio to e^(h/ε).

That is the whole of the mechanism, and it says immediately what the scheme is good for: it is exact because the answer is a geometric sequence and the scheme was built to reproduce one. Nothing about it is a general repair for convection.

Writing the quadratic out is what turns that from an observation into a list of conditions. Applying a three-point operator with diffusion ε̃ to a geometric sequence of ratio q multiplies it by

(−ε̃/h² + 1/2h)·q  +  2ε̃/h²  +  (−ε̃/h² − 1/2h)/q

and the scheme is exact when q = r is a root of that. So the derivation needs a constant velocity, so that r is one number rather than one per cell; a constant ε, for the same reason; and a right-hand side that leaves the solution a combination of 1 and rⁱ. Change any of the three and the tuned scheme is a first-order method with an unusually large constant — which is what the manufactured problem below measures, and it is the only thing this scheme can be when its assumption is not met.

Both older schemes are limits of it #

So the field’s two existing schemes are not two ideas with a third one beside them. They are the two ends of one curve, and the practical statement is that upwinding over-corrects at every finite Péclet number — by a factor of 1/ξ, which is 1.5 at Pe = 3 and 7.8 at Pe = 0.39.

That is worth stating as a criticism of upwinding rather than as praise of the tuned scheme, and it is the criticism the convection field’s own measurement already implied: upwinding’s answer is 0.388 from the exact solution of the problem it solves and 71.1 from the exact solution of the problem that was posed. It is solving a well-behaved problem that is not the one asked about, and ξ says exactly how much of that substitution was necessary.

And it cannot oscillate #

The tuned operator’s upper off-diagonal is −ε̃/h² + 1/2h, which is nonpositive exactly when ξ ≥ 1 − 1/Pe — and coth(Pe) ≥ 1 makes that true at every Pe. So the matrix is an M-matrix at every grid and every diffusion, and the solution obeys the maximum principle by inequality rather than by luck.

The last row is where central differencing crosses back under Pe = 1 and stops oscillating, which is the threshold that essay bisected to 1.0000000000000002 arriving here as a column of zeros.

The column of zeros is a consequence rather than a coincidence, and it is worth writing the chain out because it says how much diffusion the scheme is buying that property with:

ε̃  ≥  h/2       ⟺       1 + ξ·Pe  ≥  Pe       ⟺       ξ  ≥  1 − 1/Pe

Central differencing satisfies the same inequality exactly when Pe ≤ 1. The tuned scheme satisfies it at every Pe, with room to spare, because coth(Pe) ≥ 1 for every positive Pe. So the two schemes differ by whether a single condition holds at one Péclet number or at all of them — and the tuned scheme buys the second with the smallest quantity of diffusion that does it, which is what the comparison against upwinding’s h/2 in the previous table measures.

Exact at the nodes is not exact #

The nodal values are the true solution to rounding. Between them, the computed function is whatever is drawn through those points — a straight line, in every plot and in every interpolation a code would use for output.

Nearly half the solution’s whole range, in the cell nearest the boundary — because the true solution climbs from near zero to one inside that single cell, and a straight line across it is as wrong as a straight line can be.

That is not a criticism of the scheme; it is what “nodally exact” means, and the phrase is in every account of the method precisely because the distinction is real. What the scheme guarantees is the value at the points it was asked about. The layer is still unresolved, and if what is wanted is the flux at the wall, or the position of the layer, or the solution anywhere between two grid points, the tuned scheme has not helped at all.

The finding: the exactness belongs to the problem #

Same operator, same ε, same grid, same tuned ε̃. Different right-hand side: a manufactured smooth solution u(x) = sin(πx), with whatever f that requires.

The ratio grows with refinement, because central differencing is second order here — 4.00 and 4.00 to two decimals, which is asserted — and the tuned scheme is not: 3.03 and 3.61, between first and second order, because the diffusion it adds is only O(h²) in the limit Pe → 0 and these grids are not there.

The added diffusion is an error wherever there is no layer for it to stabilise, and nothing in the formula knows whether there is one. ξ is a function of Pe alone. It does not look at the right-hand side, and the right-hand side is what decides whether the solution is a geometric sequence.

Which is a finding this site has made before, in another field #

A multigrid smoother built for anisotropy returns a convergence factor of 0.037 on the problem it was aimed at and 0.967 on the same problem turned through a right angle. The name that phase gave the shape was where a method’s competence was put into it by hand — the repair encodes not “cope with anisotropy” but “cope with anisotropy along y”, and the direction was supplied by whoever typed the string.

This is the same shape with the hand-supplied information being the form of the solution rather than a direction. ξ encodes “the answer is a geometric sequence with ratio e^(h/ε)”. On a problem where that is true, the scheme is exact. On a problem where it is not, the scheme is a first-order method competing against a second-order one, and losing by a factor that grows as the grid is refined.

What a real code does with this #

The one-dimensional streamline diffusion method is this scheme — the streamline-upwind Petrov–Galerkin stabilisation reduces to exactly ε̃ = ε(1 + ξPe) on this problem — and in more than one dimension it is where the method’s design decisions actually live. Two of them are visible from here.

The direction. Streamline diffusion adds its diffusion along the flow rather than isotropically, which is what stops it from smearing the solution across the streamlines. In one dimension there is only one direction and the distinction is invisible; in two it is the whole method, and it is the same distinction the anisotropy field spends four essays on.

And the parameter. Every practical code uses ξ, or a cheaper approximation to it, computed per element from a local Péclet number. What this essay says about that is not that it is wrong but that its optimality is inherited from a one-dimensional constant-coefficient problem with no source term — and that the measurement above, of the same formula on a smooth problem, is what the inheritance costs when the local situation is not the one it was derived from.

The formula cancels, and it costs nothing #

ξ = coth(Pe) − 1/Pe is a subtraction, and coth(Pe) → 1/Pe as Pe → 0. So at small Péclet numbers the formula computes a quantity of size Pe/3 as the difference of two quantities of size 1/Pe, which is exactly the arithmetic this collection’s first field exists to be suspicious of.

It is as bad as that description makes it sound. Against the series ξ = Pe/3 − Pe³/45 + 2Pe⁵/945, the naive formula in double precision is wrong in its eighth significant digit at Pe = 10⁻⁴, in its fifth at 10⁻⁶, in its second at 10⁻⁷, and at Pe = 10⁻⁸ it returns exactly zero for a quantity whose value is 3.33·10⁻⁹. Every digit is gone, and nothing about the expression warns anybody: it is two library calls and a subtraction, and both library calls are correctly rounded.

And it does not matter at all, for a reason that is exact rather than lucky.

ξ never enters the scheme alone. It enters as ε̃ = ε(1 + ξ·Pe), so the 1/Pe that amplified the rounding error is multiplied straight back out by the Pe standing next to it. The absolute error in ξ is about u/Pe; the absolute error in ε̃ is that times ε·Pe, which is u·ε — the size of a single rounding of ε itself. Measured at ε = 10⁻³, the naive and series forms of ε̃ agree to within 3·10⁻²⁰ at Pe = 10⁻⁴, 10⁻⁶ and 10⁻⁸, against u·ε = 1.1·10⁻¹⁹. The cancellation is total, and it is invisible one line downstream.

That is worth stating carefully, because both halves are true and the conclusion is neither of them. A catastrophic cancellation is a statement about a quantity, not about a computation; it costs something only when the damaged quantity is then used in a way that preserves its relative error. Used additively, against something of the same size as the amplification, the damage cancels with the amplification that caused it. The site’s habit of printing the residual rather than the digits is the same distinction from the other side: what a number is worth depends on what is done with it next.

There is still a reason to use the series, and it is not accuracy. At Pe below about 10⁻⁸ the naive form returns a hard zero, so a code that branches on xi > 0 — to skip the stabilisation, say, or to assert that it is positive — takes a different branch than one that computes 3.33·10⁻⁹. The quantity is negligible either way; the control flow is not.

Two schemes separated by ε #

The other end of the range gives an equally sharp statement, and this one is about how fine a knob the stabilisation parameter is.

coth(Pe) = 1 + 2/(e²ᴾᵉ − 1), so ξ = 1 − 1/Pe + 2/(e²ᴾᵉ − 1). Past about Pe = 20 that last term is below the unit roundoff and coth returns exactly 1, so the tuned scheme computes ξ = 1 − 1/Pe exactly and its added diffusion is ξ·h/2 = h/2 − ε. Upwinding adds h/2. The two schemes differ by exactly ε, to within 2εPe/(e²ᴾᵉ − 1), and that correction is 6·10⁻⁵ at the Péclet number of 3.125 in the first table.

So at ε = 0.005 and n = 31 the exact scheme’s total diffusion is 0.0157 and upwinding’s is 0.0206. Those are the same number to one significant figure, they are 31 per cent apart, and one of them gives an answer correct to 2.4·10⁻¹⁷ at every node while the other is 0.136 away in the maximum norm.

That is the argument against the way stabilisation is usually described. Add a little artificial diffusion to damp the oscillation sounds like a coarse knob, tolerant of the amount, and it reads that way because the visible failure it repairs — the wiggles — is repaired by any sufficient amount. It is not a coarse knob for accuracy. The interval between the amount that is exact and the amount an ordinary one-sided difference happens to supply is one ε wide, and the answers either side of it differ by a sixth of the solution’s whole range. Getting the oscillation out is easy and says nothing about whether the answer is right, which is the same lesson the convection essay reached from the other direction when it found a well-behaved answer to the wrong problem.

What a code should report, given all of that #

Three of the measurements above point at the same practical conclusion, and it is not the one the formula’s elegance suggests.

The scheme is nodally exact on the problem it was derived for, first-order-with-a-large-constant on a problem with a source term, and its added diffusion is within ε of upwinding’s at every Péclet number above about three. So the quantity that decides whether it is helping is whether the solution is a geometric sequence, and nothing in the code knows that. A residual is no help: the tuned scheme solves its own system to 4.6·10⁻¹⁸ whether or not that system is the right one, which is the distinction between solving accurately and answering correctly that this collection returns to more often than any other.

What is available is the comparison. Running the same problem at ε̃ and at ε and reporting the difference between the two solutions costs one extra solve and measures exactly the thing at issue: how much of the answer is the stabilisation rather than the equation. On the layer problem that difference is the 0.517 error central differencing carries, and the tuned answer is the correct one. On the manufactured smooth problem it is 6.7·10⁻² against central differencing’s 1.6·10⁻³, and the untuned answer is the correct one. The number is the same measurement in both cases and only its interpretation changes, which is what makes it worth printing rather than deciding in advance.

That is a cheaper version of the residual-based schemes named below, and a much weaker one: it says how much the stabilisation moved the answer without saying whether it should have. It is still better than the alternative on offer, which is to trust a parameter derived for a problem the code is not necessarily solving — the failure that the essay on a smoother that cannot see a direction describes in a different field and the same shape.

What is left #

Two dimensions, where the flow has a direction, the stabilisation is anisotropic by design, and none of the closed forms above survive. That is where the subject actually lives and it is a phase’s work rather than a section’s.

Higher-order stabilisations — the discontinuous Galerkin family, and the residual-based schemes that add diffusion proportional to how badly the equation is being violated locally. The second of those is the direct answer to this essay’s finding: a scheme that measures whether there is a layer rather than assuming one.

And a variable coefficient. The exactness rests on ε and the velocity being constant, so that r is one number rather than a different ratio in every cell. A problem with a varying velocity has a different geometric ratio per cell and the scheme is no longer exact anywhere — which is measurable, and is not measured here.

A parameter with a floor rather than a direction #

The tuned parameter here is exact at one place and harmful next door, which is the shape of a knob with a best value rather than a best direction. The clearest measured instance of that shape elsewhere is the number of squarings in a matrix exponential.