The inner product the mesh already computed

A certificate written in coordinates took a stability guarantee that costs nothing — every pole of a Galerkin reduced model lies inside the numerical range of the operator, so a range in the left half plane certifies every model — and showed that it belongs to the coordinates the state is written in. Rescale the state by a diagonal matrix and the system is the same system, with the same eigenvalues and the same transfer function, but its numerical range moves; at Péclet 40 a rescaling with condition number 260 pushed the range’s edge to +6.6 and ten placements of nineteen returned unstable reduced models.

Its repair for a code that meets this was the energy inner product. “For a discretised physical operator those are usually the coordinates of its energy norm, which a code knows before it builds anything. The failure measured here is what happens when the state is stored in some other coordinates and the projection is done in them anyway.” And it noted why its own original coordinates had been safe: “On a uniform grid that norm is the Euclidean one scaled by the mesh width.”

That leaves a gap between the two sentences, and it is the gap every real code sits in. The rescalings that broke the certificate were constructed, powers of a ratio chosen to skew the operator. Nobody stores a state that way. What codes store is nodal values, and on a mesh that is not uniform the nodal vector is not the energy norm’s coordinates either.

A mesh graded towards the layer #

The operator is the same centred-difference discretisation of $\varepsilon u'' - u'$ on the unit interval, thirty interior nodes, with $\varepsilon$ the reciprocal of the Péclet number. At high Péclet number its solutions have a boundary layer of width about $\varepsilon$ at the outflow end, and resolving that layer is the whole reason to grade a mesh: small cells at the outflow, large ones upstream, the cells growing geometrically so that the largest is a stated multiple of the smallest.

A uniform mesh of thirty-one cells has no node inside a layer one hundredth wide. Graded a hundred to one it has seven, which is the point. On a nonuniform grid the centred stencils change shape: the second difference at node i divides by both neighbouring cell widths and their sum, and the first difference divides by the sum. Neither is symmetric or skew as a matrix any more, and the input and output vectors are left as the uniform system had them, from the same seed, so that nothing about the comparison moves except the mesh.

The nodal range crosses at a ratio of ten #

In nodal coordinates the range’s edge is −0.10 on the uniform mesh at Péclet 100, the diffusion term’s own margin. Grade the cells three to one towards the outflow and it is +0.48; ten to one, +1.38; a thousand to one, +8.94. At Péclet 40 it crosses at a ratio of three and reaches +4.18; at Péclet 10 it crosses at ten and reaches +1.19. The certificate is gone on every mesh anyone would use to resolve the layer.

Grade the cells the other way, small at the inflow where nothing happens, and the edge moves the other way: at Péclet 100 it falls to between −0.84 and −0.59, further left than on the uniform mesh. The grading is not harmful as such. It is harmful in the direction the physics asks for.

The third line is the projection in a different inner product, and it is the reason this is worth measuring. It sits at −0.098 at every grading at Péclet 100, −0.24 at Péclet 40 and −0.97 to −0.99 at Péclet 10, which are the diffusion term’s own eigenvalues on the uniform mesh to two figures. Weighted by cell size, a graded mesh has the numerical range of the continuous problem.

Why a diagonal restores it #

Multiply row i of each stencil by the width of the dual cell around node i, half the sum of its two neighbours. The second difference becomes the stiffness stencil — one over the left width, minus the sum of the two reciprocals, one over the right width — which is symmetric and negative definite whatever the widths are. The first difference becomes $(u_{i+1} - u_{i-1})/2$, which is exactly skew-symmetric. So with $M$ the diagonal of dual-cell widths, $MA$ is a symmetric negative definite matrix plus a skew one, and in the inner product $\langle x, y\rangle = x^{\mathsf T}My$ the operator is dissipative on every mesh: the convection conserves energy and the diffusion removes it, as in the continuous equation.

That diagonal is not an extra computation. It is the lumped mass matrix of a linear finite-element discretisation on the same nodes, and a finite-difference code on a nonuniform grid computes the same widths to form its stencils. A Galerkin projection in that inner product replaces $V^{\mathsf T}AV$ by $(V^{\mathsf T}MV)^{-1}V^{\mathsf T}MAV$ on the same basis $V$: the same rational Krylov space, built from the same shifted solves, projected differently. It costs a diagonal scaling of the basis and an r-by-r solve.

Where nodal models actually fail #

A range in the right half plane removes the certificate. It does not make a model fail, and a certificate written in coordinates found a whole direction of rescaling where the edge reached +6.6 and no model failed. So the count of unstable models has to be measured separately.

At Péclet 10 none do, at any grading, although the edge is past +1 from a ratio of three hundred. At Péclet 40 two models fail at a ratio of three hundred and four at a thousand. At Péclet 100 four fail at a ratio of thirty, five at a hundred and seven at three hundred and a thousand. Nearly all the failures are at order two — six of the seven at the worst mesh — with one at order six and none at orders one or three. The weighted projection returns no unstable model on any of the forty-two meshes of the sweep, at any order.

So the certificate is lost well before any model fails: at Péclet 100 the edge is past zero from a ratio of three, and the first failure is at thirty, a factor of ten further along the grading. That gap is also why a code can run for a long time on a graded mesh without meeting the problem. Coarse grading at moderate Péclet number removes the guarantee and leaves the models stable, and the guarantee is not something a code checks.

One basis, three projections #

At the worst mesh measured, the nodal order-two model is unstable at the six lowest placements, from s = 0.10 to 0.68, with its rightmost pole climbing to +1.23. On exactly the same basis, projected with the cell-size weight, the rightmost pole never rises above −1.04. The two-sided model, which matches slopes as well as values and has no certificate of any kind, is unstable at the four highest placements, from s = 31.6 up, reaching +4.7 — the failure a model that cannot be run found on the uniform mesh, still there, at the other end of the axis.

That figure is the argument in one picture. The rational Krylov space is a property of the shifts and the operator, and the nodal and weighted lines were computed from the same one. Everything that differs between them is a choice of inner product in the last step, a choice a code makes whether or not it knows it has.

What the weight costs in accuracy #

A certificate bought with worse models would be a poor trade, and half the conditions and a certificate had already found the one-sided construction giving up the slope conditions for its guarantee. The weighted projection does not change which conditions are matched — a Galerkin projection onto a rational Krylov space interpolates at its shifts in any inner product — so the question is only whether the error away from the shifts moves.

It moves in the weight’s favour. At Péclet 100 on a mesh graded a hundred to one, the order-six weighted model’s relative H∞ error is 0.19 at s = 0.1 against the nodal model’s 0.45, and 0.18 against 0.37 at s = 1. At the high placements, from s = 10, the two are within a few per cent of each other and both near one: order six on shifts that high does not capture this system’s low-frequency behaviour in either inner product. At no placement is the weighted model worse by more than the width of the rounding in an H∞ estimate.

That these errors are large at all is a fact about this system and this order, not about the weight: a transfer function dominated by a sharp layer at the outflow needs more than six directions. What the figure settles is that the certificate is not bought with accuracy. On the placements where the reduced model is worth having, the inner product the physics supplies is also the one that approximates better.

What the nodal edge is measuring #

A numerical range in the right half plane for a stable operator is a statement about transient growth, and that is what makes it more than a technicality. In the coordinates where the edge is positive, some initial states grow before they decay: the fastest rate at which the Euclidean norm of a state can grow at an instant is the edge times the norm. A spectral radius that grows first measured the discrete version of the same thing, powers of a matrix whose spectral radius was below one rising by a factor of twenty thousand before they turned over, and the eigenvalues that are not there drew the resolvent regions such a matrix carries around its spectrum.

On the graded mesh that growth is not physical. The continuous equation’s energy never grows, and in the weighted inner product the discrete one’s does not either. What grows is the Euclidean norm of the nodal vector, which on a graded mesh weighs the many small cells at the outflow as heavily as the few large ones upstream, so that a state concentrated in the layer counts for more than its energy. The nodal range’s edge is the rate at which the solution can move mass from large cells into small ones, measured in a norm that counts nodes rather than mass. A Galerkin projection that uses that norm is projecting a question about nodes, and its models inherit the transient growth of the question.

That also explains the direction. With the small cells at the inflow, convection carries mass from small cells into large ones, the nodal norm shrinks faster than the energy, and the edge moves left. With them at the outflow it carries mass the other way. The grading the physics asks for is the grading under which counting nodes overstates the solution most.

Against the other two weights #

The rescaling essay measured two other ways to restore the certificate. The symmetrising scaling exists only for an operator similar to a symmetric one by a diagonal, and its conditioning grows geometrically with the grid and the Péclet number, so that on a fine grid at high convection it cannot be represented in the arithmetic it is meant to certify. The Lyapunov solution is a dense matrix equation, the cost that interpolatory reduction was adopted to avoid, and the solver there declined to produce it above about Péclet 50.

The cell-size weight has neither limit. It is diagonal, its entries are cell widths, which differ by the grading ratio and not by anything geometric in the grid size, and it certified every model measured at Péclet 100, where the Lyapunov weight was not available at all. Its limitation is the one the others do not have: it certifies because this operator is a symmetric diffusion plus a skew convection in the energy inner product, and it would not certify an operator whose dissipative structure lives in some other norm.

What a code that reduces a discretised model should do #

Project in the inner product the discretisation defines, not in the Euclidean inner product of whatever vector the state is stored as. For a finite-element model that is the mass matrix, which is what a Galerkin reduction of such a model usually uses anyway. For a finite-difference or finite-volume model on a nonuniform grid it is the diagonal of dual-cell widths or cell volumes, which were computed to form the stencils and which it is easy to forget in the reduction.

And where the physics’ inner product is not known, compute the numerical range’s edge in the inner product actually used — a symmetric eigenvalue problem of the full size, cheap next to the reduction itself — before trusting the certificate. On the nodal graded meshes here that check would have said “no certificate” at a ratio of three, long before a model failed, which is exactly the warning a code needs.

The two remedies are not alternatives. The check tells a code whether its present inner product certifies anything; the weight gives it one that does. A code that applies the weight and then checks the weighted edge is checking a quantity it knows in advance to be the diffusion term’s eigenvalue, which is a useful test of its own implementation and nothing more. A code that checks without weighting learns, on these meshes, that it has no guarantee — and has to decide what to do about models that are, most of the time, stable anyway. Of the forty-two meshes measured — seven gradings, two directions, three Péclet numbers — the nodal projection returned an unstable model on six and stable models with no certificate on eleven more. The weight turns all seventeen into certified ones, at no cost in accuracy that the H∞ estimates can see.

What this does not settle #

One dimension, one discretisation, thirty nodes, one grading law. A two- or three-dimensional mesh graded towards a wall has a dual-cell weight that is still diagonal and a stiffness matrix that is still symmetric under it, so the argument carries; the size of the nodal edge’s excursion, which depends on how the stencils’ asymmetry adds up over neighbours, is not measured.

The failure count is a count at nineteen placements and four orders, and why order two fails so much more often than orders one, three and six is not explained. It matches the rescaling experiment, where failures also concentrated at particular orders, and it is a reminder that a lost certificate permits failure without predicting where.

The account of the nodal edge as mass counted by nodes — convection carrying the solution from large cells into small ones, which the nodal norm overweights — is an interpretation that fits the direction and the growth with Péclet number. It is not a measurement: the transient growth itself, the peak of the nodal norm over time for the worst initial state, is not computed here, and it is the quantity that would test the reading.

The H∞ errors are estimates from a frequency grid of 80 points refined around the maximum, and differences of a few per cent between the two projections at the high placements are inside that estimate’s accuracy.

Still open: upwinding, and a certificate for the two-sided model #

An upwind discretisation. Every operator here is centred, and centred differences at a cell Péclet number above two oscillate. A different equation on every grid showed that upwinding is the exact discretisation of the same equation with diffusion ε + h/2, and a parameter that is also a price measured what the fitted fraction of that extra diffusion costs on a problem with no layer. A code that grades its mesh to resolve the layer usually also upwinds the convection, which adds artificial diffusion. Upwinding changes the skew part into a skew part plus a dissipative one in the weighted inner product, and whether that makes the nodal range safe again — or merely moves where it crosses — is the next measurement on the same meshes.

The two-sided model at high placements. The one-sided construction now has a certificate in the right inner product on any grading. The two-sided one never had one and still fails at the high placements, from s = 31.6 at Péclet 100 on the graded mesh. Whether a two-sided projection whose left space is chosen in the weighted inner product — a Petrov–Galerkin projection that treats $M$ as the natural pairing — inherits any of the certificate is a question the same measurement can answer.