MEng project · UTIAS · with Prof. David Zingg
A provably stable scheme needs a norm. Build it from a quadrature rule and nothing is left to choose.
A summation-by-parts operator carries a diagonal norm, and the diagonal of that norm is also the weight vector of a quadrature rule. This project builds it as a Gregory end-corrected quadrature rather than from an ansatz, and asks whether that route produces anything new. At minimal span it has exactly one free number and one condition to spend it on.
Building the norm from a quadrature
The routeIntegration by parts is the identity that makes an energy estimate work. If a discretisation reproduces it exactly, stability is a matter of algebra rather than of hope. Summation by parts is that identity, discretely: a difference operator D and a positive norm H such that H D + (H D)T = E, where E is nothing but the boundary terms.
The route below builds H from quadrature rather than from an ansatz: the diagonal of a norm that admits an energy estimate is also the weight vector of a quadrature rule, and Gregory's end corrections are the classical way to build one. The question is whether that route produces anything new.
One letter, used throughout. These operators come in a family indexed by p: the difference is accurate to order 2p in the interior and to order p in the few rows at each end, which is why the classical ones are named 2–1, 4–2 and 6–3. Everything here is p = 1, 2 and 3. The norm’s quadrature rule is exact to degree 2p − 1, and the number of corrected weights at each end is 2p: two, four and six. Three more, because the algebra below cannot be read without them. u and v are any two grid functions, which is what makes the first line an identity rather than a special case. λ is the diagonal of H written as a vector, so the eighteen numbers in the figure above are λ. And μ is the one direction the weights are free to move in: the finite-difference vector that spans the null space of the exactness conditions, so μ(2p−1) is the one belonging to order p, not a power of anything. A superscript T is a transpose throughout, and N is the number of grid points a check is run at, which is why the verifications say 12 and 8 while the figure draws 18.
what cancels
The 2–4 operator's first row of D is −24/17, 59/34, −4/17, −3/34: denominators of 17 and 34 that appear nowhere in the problem statement. They come from H00 = 17/48.
Multiply through and they are gone: the first row of Q is −1/2, 59/96, −1/12, −1/32. The 17s were never structural. They were the norm showing through the operator.
why exactly 2p − 1
The interior weights are 1, so all the information is in the r boundary weights at each end. Requiring exactness for xk, k = 0…2p−1, is 2p linear conditions on r = 2p unknowns.
They are not independent: the system has a one-dimensional null space, spanned by the finite-difference vector μ(2p−1). So the weights are pinned to a line, not a point, and pushing to degree 2p makes the system inconsistent, which is why the family stops where it does.
why one number and not two
Everything before this step has already been spent. Q + QT = E fixes the symmetric part of Q. Exactness to degree 2p − 1 is 2p conditions on 2p unknowns whose matrix drops rank by one, which is where the free scalar γ comes from.
So exactly one requirement remains unimposed, it is linear in the weights, and a linear condition on a line is a division. That is the whole of the uniqueness result: not that no other norm exists, but that at minimal span there is no second number to vary.
Uniqueness, demonstrated by breaking it
The experimentThe claim is that one value of γ works and no other does. Below, λ = σ + γ μ(3) is rebuilt at five chosen values of γ, and three separate properties are recomputed each time in exact rational arithmetic: positive-definiteness of H, the summation-by-parts residual Q + QT − E, and the degree to which the quadrature is exact. Nothing here is evaluated in floating point.
γ is not a parameter with a range of admissible settings. It has exactly one admissible value, and the useful question is what specifically fails at the others, which is a question with discrete answers.
which value of γ
γ =
| k | quadrature | exact |
|---|
Narration
01 The closure value
At γ = 1/48 the three properties hold at once. Every weight is positive, so H is a norm; the summation-by-parts residual Q + QT − E is zero in all 144 entries rather than merely small; and the quadrature is exact to degree 3, which is 2p − 1.
02 One step away
Set γ to zero and you have the raw Gregory seed, before the closure is applied. The quadrature does not notice: still exact to degree 3, still failing at 4. What breaks is the identity, in 34 entries, the largest of them 51/688. Twice the closure value fails by the same amount with every sign reversed: the same 34 entries, the same 51/688. Q is linear in H, so the residual is linear in the distance from the closure, and one step either side of it is the same size of error either side of zero.
03 Far enough to lose the norm
At γ = 25/48, twenty-four steps out, the residual is twenty-four times larger at 153/86, and λ₀ and λ₂ have changed sign. H is no longer a norm, so there is no energy estimate left to make. The same distance the other way puts λ₁ negative and returns the same 153/86. The two failures are not one failure arriving twice: the identity goes immediately and in proportion to the distance, and positivity goes only once a weight crosses zero.
What the experiment shows
The quadrature is exact to degree 3 = 2p − 1 at every value of γ, and fails at degree 4 at every value of γ. Accuracy does not select the norm: the whole one-parameter family integrates equally well, which is exactly what it means for μ to span the null space of the moment conditions.
What selects it is compatibility with D, which is what the run above walks. The summation-by-parts residual is zero at one value of γ and at no other, and away from it the residual grows in proportion to the distance, so the closure is not the best fit among candidates. It is the only member of the family that satisfies the identity at all.
The same answer, from six different seeds
Robustness of the route| p, the boundary order | seed | σ differs from λ by | closure γ | resulting norm |
|---|
The route works, and it recovers what was already there
The verdictThe route works, and it produces the operators that were already there. At minimal span the quadrature view leaves exactly one scalar free and the closure consumes it, so the diagonal norm is forced, which means this construction recovers the classical Kreiss–Scherer and Strand families rather than extending them. A uniqueness theorem is a statement about what cannot be built, and that is the result. New operators would have to come from widening the stencil past minimal span, where the null space is no longer one-dimensional and there is something left to choose.
Scope
- Diagonal norms, first derivative, minimal span r = q = 2p. Nothing here applies to block norms or to wider stencils.
- The operators shipped in SBPLab are the published Kreiss–Scherer (1974) / Strand (1994) family, used as such. The Gregory construction is in the MEng report and is not implemented in SBPLab; the two should not be read as one body of code.
- Uniqueness is subject to H ≻ 0, which is checked, not assumed.
- The report prints the p = 3 closure constant as 1399687/30. Two independent routes give the same value, one from its own γ and its own weights and one backwards from its own published λ: 1399687/180. Nothing downstream depends on it, and the constant above is computed rather than quoted.
What I would do differently
Push to p ≥ 3 earlier. Everything at minimal span is forced, and that is the result, so the interesting engineering only starts where the closure stops being unique and the residual parameters have to be chosen against competing criteria. I spent most of the project in the region where the algebra decides, and the region where a person decides is the one worth a longer look.
The project closed with the MEng in 2025. No manuscript followed and none is planned, but the questions carried on: they continue in SBPLab, a private C++ bench for SBP-SAT operators and canonical problems with exact answers, the same way the CFD coursework carried on as CFDLab.