Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

The diffusion that makes the answer exact

Upwinding adds h/2 of artificial diffusion. Central differencing adds none. Add ε·ξ·Pe with ξ = coth(Pe) − 1/Pe and the computed solution is the exact one at every grid point, to 2.4·10⁻¹⁷ — at every Péclet number, on the problem it was derived from and on no other.
17 min read 9 figures Exact ground truthWhere the algebra stops being the arithmeticMeasured, not asserted

Worth reading first: The stencil that is not symmetric.

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.

Three schemes at ε = 0.005, cell Péclet number 3.13The computed solution of each scheme against position, with the exact solution drawn as a dashed line. The tuned scheme's values lie on it — the largest nodal difference is 2.4·10⁻¹⁷. Upwinding is smooth and 0.136 away at its worst. Central differencing oscillates and leaves the interval [0, 1] at 16 of the 31 points.00.250.50.75100.51xuexact & tunedupwindcentralone of these is exacttuned, worst nodal error2.4·10⁻¹⁷upwind, worst nodal error0.14central, points outside [0, 1]16exact at every nodeand only at the nodes
Fig. 1 Three schemes at ε = 0.005 on a 31-point grid, cell Péclet number 3.125, with the exact solution dashed. The tuned scheme’s values lie on it — the largest nodal difference is 2.4·10⁻¹⁷. Upwinding is smooth and 0.136 away at its worst; central differencing oscillates and leaves [0, 1] at sixteen points. Drag ε and watch which curve stays on the dashed one.

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.

n ε Pe ξ tuned upwind central
31 0.005 3.125 0.684 2.4·10⁻¹⁷ 0.136 0.517
31 0.05 0.313 0.103 5.6·10⁻¹⁷ 0.092 0.012
63 0.002 3.906 0.745 2.7·10⁻¹⁷ 0.113 0.593
15 0.1 0.313 0.103 1.1·10⁻¹⁶ 0.092 0.012

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.

Both older schemes are limits of it

The optimal upwinding parameter, and the two schemes it containsξ = coth(Pe) − 1/Pe against the cell Péclet number on a logarithmic horizontal axis. It runs from zero to one: at small Pe it behaves like Pe/3, so the added diffusion vanishes faster than the grid and the scheme becomes central differencing; at large Pe it approaches one — 0.9667 at Pe = 30 — where the added diffusion is h/2 and the scheme becomes upwinding.10⁻²10⁻¹110¹00.250.50.751cell Péclet numberξupwinding: ξ = 1Pe = 1ξ(Pe)Pe/3one curve, two schemesξ at Pe = 0.10.033ξ at Pe = 10.31ξ at Pe = 300.97central differencing at one endupwinding at the other
Fig. 2 ξ against the cell Péclet number. At small Pe it behaves like Pe/3 — checked against ⅓ to 2% at Pe = 10⁻², 10⁻³ — so the added diffusion vanishes faster than the grid and the scheme becomes central differencing. At large Pe it approaches one like 1/Pe, so the added diffusion approaches h/2 and the scheme becomes upwinding.
n Pe ξ tuned adds upwind adds
31 3.125 0.684 1.07·10⁻² 1.56·10⁻²
63 1.563 0.452 3.53·10⁻³ 7.81·10⁻³
127 0.781 0.250 9.78·10⁻⁴ 3.91·10⁻³
255 0.391 0.129 2.52·10⁻⁴ 1.95·10⁻³

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.

n Pe tuned: points outside [0, 1] central: points outside
7 12.5 0 4
15 6.25 0 8
31 3.125 0 16
63 1.563 0 11
127 0.781 0 0

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.

−εu″ + u′ = 0 on 31 points, ε = 0.005, cell Péclet 3.125Three solutions of a boundary-layer problem on the same grid. The exact one rises monotonically from 0 to 1 with a layer of width ε at the right-hand end. The central-difference answer alternates at every point and leaves [0, 1] at 16 of the 31 of them, by as much as 0.5152. The upwind answer is monotone at every Pe.00.250.50.75100.250.50.751xuexactcentral differencesupwindthe oscillation is exact‖Ax − b‖/‖b‖ for the central answer4.6·10⁻¹⁸values outside [0, 1]16worst excursion0.52the dashed lines are 0 and 1, which the equation guaranteesno solver was involved
Fig. 3 The threshold itself, from the field’s own essay: the off-diagonal −ε/h² + 1/2h changes sign at exactly Pe = 1, and that is where the matrix stops being an M-matrix and the answer starts leaving the interval. The tuned scheme’s off-diagonal cannot change sign at any Pe, which is the one-line reason its column above is zeros.

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.

n ε nodal error midpoint error where
31 0.005 2.4·10⁻¹⁷ 0.457 x = 0.984
63 0.002 2.7·10⁻¹⁷ 0.480 x = 0.992

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 same three schemes on a problem with no layer in itWorst nodal error against the grid size, both axes logarithmic, for a manufactured smooth solution on the identical operator at ε = 0.005. Adding no diffusion gives 0.0016, 4·10⁻⁴, 10·10⁻⁵, falling by four at each refinement. The tuned diffusion gives 0.067, 0.022, 0.0061 — 42 times worse at the coarsest grid, and falling more slowly.10²10⁻⁴10⁻³10⁻²10⁻¹grid points nworst nodal errorupwindtunedcentralthe same tuning, another problemtuned ÷ central at n = 3142tuned ÷ central at n = 12761central's error at the finest grid10·10⁻⁵exact on the problem it was derived fromand harmful on the one beside it
Fig. 4 The three schemes on a problem with no layer in it. Adding no diffusion gives 1.60·10⁻³, 3.99·10⁻⁴ and 9.96·10⁻⁵ as the grid is refined — falling by exactly 4.00 each time, which is second order. The tuned diffusion gives 6.69·10⁻², 2.21·10⁻² and 6.12·10⁻³: forty-two times worse at the coarsest grid, and falling more slowly.
n tuned central tuned ÷ central
31 6.69·10⁻² 1.60·10⁻³ 41.9
63 2.21·10⁻² 3.99·10⁻⁴ 55.4
127 6.12·10⁻³ 9.96·10⁻⁵ 61.4

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

The depth phase measured a multigrid smoother built for anisotropy: it 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.

The same anisotropy along each axis, ε = 0.001Six horizontal bars of convergence factor. The pairs aimed at the strong direction are short; the same methods on the same problem turned sideways are nearly full length.point · strong in y0.951point · strong in x0.956y-line · strong in y0.037y-line · strong in x0.967semi-y · strong in y0.107semi-y · strong in x0.9661.00 — no convergenceone problem, seen from two sidesy-line, aimed0.037y-line, turned sideways0.97the ratio2631×31 grid, twelve V-cyclesthe direction is in the code, not in the problem
Fig. 5 The multigrid version of the same sentence, from the depth phase: two repairs that work on the problem they were aimed at and return 0.967 and 0.966 on the same problem with its axes swapped. What was put in by hand there was a direction; here it is the shape of the answer.

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 upwind answer, against the two problems it might be answeringThe upwind solution drawn against two exact solutions: the one for the ε that was asked for, whose layer is thin, and the one for ε(1 + Pe) = 0.02063, whose layer is 4.125 times thicker. It lies on the second. The two operators are identical entry for entry, to 1.3·10⁻¹⁶.00.250.50.75100.250.50.751xuexact at ε = 0.005exact at ε(1 + Pe)the upwind answeran exact answer to a different question‖upwind(ε) − central(ε(1+Pe))‖/‖·‖1.3·10⁻¹⁶distance to the problem it solves0.39distance to the problem posed71the added diffusion is h/2 = 0.01563, whatever ε isso refining removes it
Fig. 6 The identity this essay’s scheme generalises: upwinding is central differencing at ε(1 + Pe), entry for entry to 1.3·10⁻¹⁶. Replace the 1 by ξ and the scheme becomes exact — so the field’s two schemes and this one differ in a single coefficient.
The damping of each Fourier mode, ε = 0.005, ω = 0.6667The modulus of one weighted Jacobi sweep's effect on each grid mode, against the frequency. On the operator with no convection in it the effect is a real number and the oscillatory half is damped to ⅓ at ω = 2/3. With convection the symbol is complex, its imaginary part does not depend on ω, and the worst damping over the oscillatory half is 1.0937.00.7853981.57082.356193.1415900.250.50.751frequency θdamping |g(θ)|the oscillatory half →⅓ — the symmetric optimumno convectionwith convectionits imaginary parta modulus, not a valuesmoothing factor at this ω1.1best over every ω0.84the symmetric operator's, at ω = 2/30.33the imaginary part does not depend on ωso no ω removes it
Fig. 7 What the same operator does to a multigrid smoother: the Fourier symbol goes complex and no relaxation reaches the symmetric operator’s smoothing factor. The tuned discretisation above changes the operator the smoother would see, which is a connection this field has not measured.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 8 The site’s other exactly known answers, for comparison. Those are formulas for an answer; ξ is a formula for a scheme, and its exactness is of the discretisation rather than of the solve.
−εu″ + u′ = 0 on 31 points, ε = 0.02, cell Péclet 0.7813Three solutions of a boundary-layer problem on the same grid. The exact one rises monotonically from 0 to 1 with a layer of width ε at the right-hand end. Below Pe = 1 the central-difference answer is monotone and inside [0, 1] like the exact one. The upwind answer is monotone at every Pe.00.250.50.75100.250.50.751xuexactcentral differencesupwindthe oscillation is exact‖Ax − b‖/‖b‖ for the central answer2.5·10⁻¹⁷values outside [0, 1]0worst excursion0the dashed lines are 0 and 1, which the equation guaranteesno solver was involved
Fig. 9 The threshold at a milder diffusion, where central differencing is still an M-matrix. The tuned scheme is one at every ε and every grid, which is what makes its column of excursion counts all zeros.

The derivation, in the form that says what it assumes

The exactness reads as luck until the two lines behind it are written down, and the two lines are also the list of what has to be true for it to hold.

The exact solution sampled at equally spaced points is a geometric sequence plus a constant: uᵢ = (rⁱ − 1)/(rⁿ⁺¹ − 1) with r = e^{h/ε}. Apply a three-point difference operator with diffusion ε̃ to a geometric sequence with ratio q and the result is that sequence multiplied by

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

which is zero exactly when q is a root of a quadratic whose coefficients are the stencil’s. So the discrete equation is satisfied exactly if the operator’s quadratic has a root at r — and ξ is the diffusion that puts it there. That is the whole mechanism.

What the derivation needs is now visible: a constant coefficient, 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 scheme is a first-order method with an unusually large constant, which is what the manufactured problem measures.

Why it cannot oscillate, as an inequality

The M-matrix property is checked as a measurement above and it holds for a reason that needs no measurement. The upper off-diagonal is −ε̃/h² + 1/2h, which is nonpositive exactly when

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

and ξ = coth(Pe) − 1/Pe with coth(Pe) ≥ 1 for every positive Pe. So the inequality holds at every Péclet number with room to spare, and the column of zeros in the excursion table is a consequence rather than a coincidence.

That is worth having beside the central scheme’s threshold. Central differencing satisfies the same inequality exactly when Pe ≤ 1, which is the threshold that field bisected to 1.0000000000000002 — so the two schemes differ by whether a condition holds at one Péclet number or at all of them, and the tuned scheme buys that with the smallest amount of diffusion that does it.

And where a real code puts this

The one-dimensional streamline diffusion method is this scheme: the SUPG stabilisation of the one-dimensional convection–diffusion equation reduces to exactly ε̃ = ε(1 + ξPe), which is why the formula appears in every finite-element account of stabilisation as well as in every finite-difference one.

In more than one dimension the method’s design decisions become visible, and two of them are worth naming from here because both are invisible in one dimension. The stabilisation is added along the streamline rather than isotropically, which is what stops it smearing the solution across the flow — a distinction that does not exist when there is one direction. And the parameter is computed per element from a local Péclet number, which inherits its optimality from the constant-coefficient one-dimensional problem above and therefore inherits the limitation this essay measures.

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.

Where this goes next

00.250.50.75100.250.50.751xuexactcentral differencesupwindthe oscillation is exact‖Ax − b‖/‖b‖ for the central answer4.6·10⁻¹⁸values outside [0, 1]16worst excursion0.52the dashed lines are 0 and 1, which the equation guaranteesno solver was involved The stencil that is not symmetric Past a cell Péclet number of exactly one — measured by bisection at 1.0000000000000002 — the central-difference solution of a convection–diffusion problem oscillates from point to point and leaves the interval the equation guarantees, at 16 of 31 grid points. It is the exact solution of its own linear system, to 4.6·10⁻¹⁸. No solver was involved. 75 aggregates over 225 unknownsthe matrix chose thismean extent along y3mean extent along x1points adopted by pass two15no coordinate enters the methodand the shape follows the coupling Aggregating what the matrix calls strong The depth phase measured every method it had on the 45°-rotated anisotropic operator — 0.784, 0.883, 0.844 — and diagnosed the failure as being in the discretisation rather than in the hierarchy. Smoothed aggregation is the standard answer to anisotropy. It returns 0.789. θx (frequency across x)θy0π/2π0π/2πthe coarse grid'sunder 0.2under 0.40under 0.60under 0.80under 0.95under 1.01damping per sweeptwo routessmoothing factor, scanned1closed form130×30 frequency cellsthe marker is the mode nothing removes A direction the smoother cannot see Give the Laplacian a strong direction and multigrid stops working — from 0.2016 a cycle to 0.9565 — with every component unchanged and the condition number identical to twelve digits. The problem did not get harder. The link between the method's two halves broke.

What links here

Computed from the collection, not written here: the essays that point at this one.

Named objects

A flat tag is an object no other essay names yet.

Artificial diffusionBoundary layerConvection diffusionDiscretisation errorM matrixManufactured solutionPeclet numberStreamline diffusionUpwind differencing