The stencil that is not symmetric

CLAIM_REGISTRY.md records this as one of four areas the depth phase explicitly declined and left for a later one: convection-dominated operators, where the stencil stops being symmetric and every Fourier argument in this phase assumes it is.

This essay is that deferral taken. The operator is one term longer than the one the whole iterative field rests on:

−ε u″ + u′ = 0      on (0, 1),   u(0) = 0,  u(1) = 1

With no first derivative it is the discrete Laplacian, whose entire spectrum this site knows in closed form and whose Jacobi, Gauss–Seidel and optimal-SOR rates it checks against closed forms. Add the u′ and every one of those closed forms stops applying, because all of them are derived from a symmetric operator with real eigenvalues and orthogonal eigenvectors, and this one has none of the three.

The oscillation belongs to the discretisation #

The claim this essay exists for, and the reason it belongs on this site rather than in a numerical methods textbook.

At ε = 0.005 on 31 points the cell Péclet number is 3.125, and the central-difference solution comes back alternating in sign from one grid point to the next — 100% of consecutive differences change sign — and leaves the interval [0, 1] at sixteen of the thirty-one points, by as much as 0.515.

The continuous problem obeys a maximum principle with no room for interpretation. Its solution rises monotonically from 0 to 1 and every value is inside that interval. The discrete answer has produced something the equation forbids.

And it is the exact solution of the system it came from: ‖Ax − b‖/‖b‖ = 4.6·10⁻¹⁸. So no amount of iteration, preconditioning, refinement or precision changes it. The badge on the figure prints that residual for exactly this reason — a reader’s first hypothesis on seeing an oscillating answer is that the solver has not converged, and the number that refutes it is the one every other figure on this site already carries.

This is the separation this collection has made since its foundation phase, moved one level further out. A small residual is not a small error separates the algorithm’s contribution from the problem’s. Here the algorithm’s contribution is at the level of rounding and the problem is not the problem either: the fault is in the discretisation, which is upstream of both and which nothing in the linear algebra can see.

The threshold is exactly one #

Not approximately one, not one within a constant. Bisecting on ε at n = 63 for the value where the oscillation appears gives a cell Péclet number of

Pe  =  1.0000000000000002

which is 1 to fifteen digits, and the assertion requires it within 10⁻⁹.

The reason is one line of the matrix. The central-difference row is

( −ε/h² − 1/2h ,  2ε/h² ,  −ε/h² + 1/2h )

and the upper off-diagonal, −ε/h² + 1/2h, changes sign at exactly h/2ε = 1. Below that both off-diagonals are negative, the matrix is an M-matrix, its inverse is entrywise non-negative, and a non-negative right-hand side gives a non-negative solution — a discrete maximum principle, inherited exactly from the continuous one. Above it the sign flips, the M-matrix property goes, and nothing guarantees anything.

The figure asserts both halves of that: the off-diagonal is negative just above the measured ε* and positive just below it. So the threshold is not fitted from the oscillation; it is read off the matrix and then confirmed by bisecting on the oscillation, which is two routes to a number applied to a threshold.

It also says which repairs work. Refining the mesh works, because Pe = h/2ε falls with h — which is the good news and is why nobody worries about this on a small problem. On a real problem the mesh that would bring Pe below one is often a mesh nobody can afford, which is why the other repair exists.

Upwinding solves a different problem, exactly #

The standard repair is to difference the first derivative one-sidedly, against the flow. It works completely: the upwind answer is monotone and inside [0, 1] at every ε, because both its off-diagonals are negative whatever h and ε are.

What it does to get there is the essay’s second measurement, and it is a backward-error statement about a discretisation:

the upwind scheme at ε   ≡   the central scheme at ε(1 + Pe)

Identically. Not approximately — ‖U − C(ε(1+Pe))‖/‖U‖ = 1.3·10⁻¹⁶, entry for entry, which is a comparison in norm only because that is how a matrix difference is reported.

So upwinding is not a better approximation of the problem posed. It is an exact discretisation of a problem nobody asked for, whose diffusion coefficient is (1 + Pe) times larger and whose boundary layer is correspondingly thicker.

The added diffusion is εPe, which is exactly h/2 — independent of ε altogether. It is a property of the mesh rather than of the equation, which is both reassuring and alarming: refining removes it, and at any fixed mesh it is whatever it is regardless of how small the physical diffusion was.

The measurement that settles it: the upwind answer is at a relative distance of 0.388 from the exact solution of the thicker-layer problem and 71.1 from the exact solution of the problem that was posed. A factor of 183.

That is the exact answer to a nearby problem — this site’s oldest idea, and the one its whole spine is built on — arriving in the one place it is almost never applied. Backward error is a statement about algorithms. Here the discretisation has a backward error: it solves a nearby equation exactly, and the nearby equation is computable.

And the multigrid field’s closed forms go with it #

A rate that is known in advance and the whole multigrid field rest on a Fourier analysis: every grid mode is an eigenvector of the sweep, the damping factor is 1 − 2ω sin²(θ/2), it is real, and the smoothing factor max(|1 − ω|, |1 − 2ω|) is minimised at ω = 2/3 where it is exactly ⅓ — a rational optimum this site found and asserted in its maturity phase.

On this operator the same calculation returns a complex number. The mode is not an eigenvector at all, and what damps is the modulus.

That is not a small change. A real damping factor can be made small by choosing ω. A modulus is bounded below by its imaginary part, and the imaginary part does not depend on ω except by an overall scaling — so no choice of ω removes it.

Scanned over forty-one values of ω from 0.05 to 1.2, the best smoothing factor achievable is 0.8425, at ω = 0.28. Against exactly ⅓ on the symmetric operator, at ω = 2/3. A smoothing factor of 0.84 is not a smoother: it is a method that removes 16% of the oscillatory error per sweep, and multigrid’s whole argument requires it to remove most of it.

The control is drawn on the same axes and is the operator with no convection in it at all, rather than one with a large ε. That distinction was a mistake the first version of this file made: setting ε = 1 leaves the u′ term in place, its contribution to the symbol is still imaginary, and the “symmetric” control measured an imaginary part of −0.005. Diffusion dominance shrinks the imaginary part and does not remove it. The symmetric operator is the one with β = 0, and on it the closed form is confirmed to 10⁻¹² before the convective case is compared against it.

The operator is genuinely far from symmetric — ‖A − Aᵀ‖/‖A‖ = 1.336 — which is the one-line explanation of why none of the field’s machinery transfers.

Both schemes are discretisations, and the repair costs an order #

Checked first, because everything above compares them and a comparison between a scheme and a bug would be meaningless. At ε = 0.5, where the layer is resolved at every size drawn, the errors fall with h at fitted orders of

central   2.011
upwind    0.988

Second order and first order, fitted rather than quoted. So the repair costs an order of accuracy, and that is the honest price beside the added diffusion: upwinding gives a monotone answer, to a different problem, at half the convergence rate.

What a reader should take from the threshold #

The practical shape is unusually clean, and it is worth separating from the mechanism.

Compute Pe before believing anything. It is h/2ε and takes no arithmetic. Below one, central differences are second-order accurate and monotone and there is no decision to make. Above one they are second-order accurate and produce answers the equation forbids.

An oscillating answer is a discretisation diagnosis, not a solver diagnosis. The first instinct is to check the solver, and the residual settles that in one line.

And the repair has a price with two parts. An order of accuracy, and a boundary layer (1 + Pe) times too thick. At Pe = 3.125 that is a layer four times too wide, which for a problem whose entire interest is the layer is not a detail.

The problem was changed to make the measurement possible #

Worth recording, because the first version measured the wrong thing convincingly.

It used the problem with f = 1 and zero boundaries, whose exact solution is a straight line minus a boundary layer. That solution has a genuine interior maximum, so its consecutive differences legitimately change sign once — and a wiggle detector built on counting sign changes then reported 3.4% oscillation on the converged, correct, non-oscillating solution at Pe = 0.03.

A measurement that is nonzero where the effect is absent cannot show a threshold, and the threshold is the whole essay.

The boundary-layer problem — f = 0, u(0) = 0, u(1) = 1 — has a monotone solution rising from 0 to 1, so a single sign change is already spurious and the detector reads exactly zero everywhere the effect is absent. It also supplies a second and better measurement for free: the maximum principle bounds every value in [0, 1], so leaving that interval is a violation that needs no judgement about how wiggly a curve looks.

Both are carried, because they measure different halves of the same defect. Near the threshold a solution can alternate without leaving [0, 1] — at n = 63, ε = 0.005 the wiggle is 0.787 and the excursion count is 11 — so a figure with only the second would show the effect arriving later than it does.

Both repairs, and which one is a repair #

Two things remove the oscillation and they are not the same kind of intervention.

Refine the mesh. Pe = h/2ε falls with h, so halving the mesh halves the Péclet number, and below one the central scheme is monotone and second-order accurate. Nothing is given up. This is the repair that is always correct and is often unaffordable: at ε = 10⁻⁶, which is an ordinary physical value, Pe < 1 needs half a million points in one dimension.

Upwind. Available at any mesh, and it changes the equation being solved. An order of accuracy and a layer (1 + Pe) times too thick.

The distinction matters because they are usually presented as two options. They are not. One resolves the problem and the other trades it for a different problem that can be resolved, and the second is only ever chosen because the first is out of reach.

Which is also why the added diffusion being exactly h/2 is the reassuring half of that measurement: it is a mesh quantity, so the artificial layer shrinks under refinement even though the scheme is first-order. The two repairs converge on each other, and the upwind scheme’s answer is wrong by an amount that vanishes — slowly, at first order, but it vanishes.

What this closes, and what it opens in the field it belongs to #

This essay sits in the iterative field rather than in a field of its own, and the placement is an argument.

Everything the field has built rests on one operator. A rate that is known in advance checks Jacobi, Gauss–Seidel and optimal SOR against closed forms derived from the discrete Laplacian’s spectrum. The error smoothing cannot reach derives the smoothing factor from the same spectrum. A rate that does not notice the size measures a V-cycle whose components were designed against it.

The depth phase already found that some of that is narrower than it looks — the Galerkin identity is one-dimensional, the anisotropy repairs encode a direction somebody typed. This essay finds the narrowest assumption of all: symmetry, which every Fourier argument in the field uses and none of them states.

What is closed is the deferral. What is opened is a question the field cannot answer with its current machinery: what a multigrid method for a non-symmetric operator looks like. The answer is known in the literature — the smoother has to be chosen for the direction of the flow, and the coarse-grid operator has to be built to respect it — and it is exactly the same shape of answer as the anisotropy repairs the depth phase measured, with the same defect: the competence is put in by hand, and the method returns a different number when the flow turns round.

That is the depth phase’s organising finding arriving in a second field, which is the sort of thing worth noticing when it happens and is not measured here.

The condition number is not where the difficulty is #

One measurement worth reporting because it repeats a finding from the depth phase in a new place.

Coarsening in one direction only found that the anisotropic problem’s condition number does not contain the anisotropy at all — both extreme eigenvalues carry the same factor, it cancels, and κ = cot²(πh/2) for every ε, measured at 103.08686891981742 across four decades with a relative spread of 2.8·10⁻¹⁶. The problem is not harder. Only the method is worse.

The same shape holds here, for a different reason and with a different quantity. The oscillation is not a conditioning failure: at ε = 0.005 on 31 points the central-difference matrix’s condition number is entirely ordinary and its solve is accurate to 4.6·10⁻¹⁸. There is nothing wrong with the linear algebra anywhere.

What is wrong is one step further out — the matrix is a faithful discrete version of a different problem from the one intended, and every quantity the linear algebra can measure is a quantity about that different problem. A residual, a condition number, a backward error: all of them are computed against the system, and the system is the thing that is wrong.

Which is the fourth level this collection has now separated. The algorithm’s fault, the problem’s fault, the method’s fault when its competence was put in by hand, and now the discretisation’s — and each of the four is invisible to the instruments that measure the ones below it.

What is left #

Two dimensions, where the flow has a direction and the mesh has two, and the interesting case is the one where they do not line up. The depth phase found that at 45° every anisotropy repair it had built fails — 0.784, 0.883 and 0.844 — because the standard nine-point stencil does not contain the anisotropy the equation has. Convection at an angle to the mesh is the same shape of problem and this essay does not reach it.

Streamline diffusion and the higher-order repairs, which add artificial diffusion along the flow only rather than isotropically, and which are the actual answer used in practice. They recover the order the plain upwind scheme gives up, at the cost of a scheme whose backward error is much harder to state.

And the solvers. Everything here is about the discretisation, and every system was solved directly. What a Krylov method does on a strongly non-symmetric operator is the GMRES essay’s subject and is not good news: this site already has a matrix with every eigenvalue on the unit circle on which GMRES makes no progress for n − 1 steps. A convergence claim for this operator from its spectrum would be exactly the claim that essay refuses.