SBP Operators

José M. Martín-García
jmm@iem.cfmac.csic.es

December 2004 (at LSU)
Updated January-March 2006 (at IEM)

Eliminate annoying spelling messages:

Off[General::spell]
Off[General::spell1]

This Mathematica notebook generates all SBP derivative operators (DOs from now on) both for first and second order, closely following the papers by Strand and Mattsson & Nordström, respectively:

"Summation by Parts for Finite Difference Approximations for d/dx", Bo Strand, Journal of Computational Physics 110, 47-67 (1994).

"Summation by parts operators for finite difference approximations of second derivatives", Ken Mattsson and Jan Nordström, Journal of Computational Physics (2004).

The structure of the notebook is as follows:
	- Section I sets the notations and analyzes the accuracy degree of a DO
	- Section II.A gives the SBP conditions for a first-order DO
	- Section II.B shows a collection of examples of different SBP DOs (those given by Strand)
	- Section III.A gives the SBP conditions for a second-order DO
	- Section III.B shows a collection of examples of different SBP DOs (those given by Mattson & Nordstrom)
	
This notebook was written during my stay at the Physics Department, LSU, in December 2004. I thank them and the Center for Computing and Technology,  LSU, for support.

For an application of this notebook see: "New, efficient, and accurate high order derivative and dissipation operators satisfying summation by parts, and applications in three-dimensional multi-block evolutions", P. Diener, E. N. Dorband, E. Schnetter, M. Tiglio, J. Sci. Comput. 32 (2007) 109-145, gr-qc/0512001.

I. Notations for derivative operators

I.1. Stencils

Let us assume a grid-function gridf on a 1d grid with points numbered from nL to nR and fixed stepsize h:

gridf[nL_, nR_] := Table[f[i h], {i, nL, nR}]

gridf[-1, 4]

Mathematica 5 gives a very general formula for a derivative stencil of maximum accuracy. There are three arguments:
	m: differential order of the derivative operator
	n: number of cells in the stencil (number of points minus 1)
	s: offset of the point of evaluation of the derivative (0 is the leftmost, n is the rightmost; centered stencils have n=2s)

UFDWeights[m_,  n_,  s_]  :=  
      CoefficientList[Normal[Series[x^s Log[x]^m, {x, 1, n}]/h^m], 
        x]

stencil[m_, n_, s_] := 
      UFDWeights[m, n, s] . gridf[-s, n - s]

Examples. We check that the accuracy degree is always n+1-m (number of points minus derivative order, independent of the offset s):

stencil[1, 4, 2]

Series[%, {h, 0, 4}]

stencil[2, 4, 0]

Series[%, {h, 0, 3}]

stencil[3, 7, 2]

Series[%, {h, 0, 5}]

In Strand we find (after correcting a few typos) the following explicit representation of the particular case m=1, n=2s:

alpha[0, _] := 0
alpha[nu_?Negative, 
        s_] := -alpha[-nu, s]
alpha[nu_?Positive, 
        s_] := -(-1)^
            nu (s!)^2/
              nu/(s + nu)!/(s - nu)!

stencilStrand[s_] := 
      Expand[Table[alpha[nu, s], {nu, -s, s}] . gridf[-s, s]/
          h]

We check it up to s=20:

And @@ 
      Table[stencil[1, 2  s, s] == stencilStrand[s], {s, 1, 20}]

I.2. Matrix derivative operators

Because DOs are linear operators, it is natural to collect the stencils for all points in a square matrix DD. We say that DD (of differential order do) has accuracy degree tau if it gives exact results on power-law grid-functions of degrees 0...tau+do-1. This condition will be imposed as a single matrix condition of the form DD.LeftMatrix == RightMatrix. We need to define a new power function to avoid indeterminacies:

power[0, _?Negative] := 0
power[0, 0] := 1
power[i_, j_] := Power[i, j]

LeftMatrix[nn, k] gives a collection of k power-law grid-functions with exponents j=0...k-1, with nn+1 points each:

LeftMatrix[nn_Integer, k_Integer] := 
      Table[power[i, j], {i, 0, nn}, {j, 0, k - 1}]

LeftMatrix[4, 5] // MatrixForm

RightMatrix gives the expected result of a derivative operator of differential order do on LeftMatrix:

RightMatrix[nn_Integer, tau_Integer, do_Integer] := 
      Table[j!/(j - do)! power[i, j - do], {i, 0, nn}, {j, 0, 
          tau + do - 1}]

RightMatrix[4, 4, 1] // MatrixForm

RightMatrix[4, 4, 2] // MatrixForm

Any derivative (of any order) will be represented by a square matrix of size nn+1:

DD[nn_Integer] := Array[d, {nn + 1, nn + 1}]

AccuracyConditions[nn_Integer, tau_Integer, od_Integer] := 
      DD[nn] . LeftMatrix[nn, tau + od] == 
        RightMatrix[nn, tau, od]

This is a system of (nn+1)(tau+od+1) equations for (nn+1)^2 unknowns. Clearly we need tau+od≤n in order to get a solution. In other words, in order to get a stencil of accuracy tau for a derivative of order od we need at least tau+od points.

This function solves the system of equations. Note that we use Solve, and not Reduce, because we assume that there are no special cases.

AccuracySolution[nn_Integer, tau_Integer, od_Integer] := 
      Print["Accuracy order is too high."] /; 
        tau > nn
AccuracySolution[nn_Integer, tau_Integer, 
        od_Integer] := Solve[AccuracyConditions[nn, tau, od], 
          Flatten@DD[nn]][[1]]

DD[3] /. AccuracySolution[3, 2, 2] // MatrixForm

Note that we do not include powers of h in the definition of the operator DD.

I.3. Checking function

This function checks the accuracy of an operator on nn+1 points, using Taylor expansions up to order taylor:

check[operator_, nn_Integer, taylor_Integer, od_Integer] := 
      Module[{gf, g, m, i}, 
gf = 
          gridf[1, nn + 1]; 
g[expr_, {i_}] := 
          expr /. f[m_.  h] :> 
              f[x + (m - i) h]; 
Take[#, 2] & /@ 
          Rationalize[
            Normal /@ Simplify@
                Map[Series[#, {h, 0, taylor}] &, 
                  MapIndexed[g, operator . gf/h^od]]]]

I.4.a. Examples: first-order systems

In the spectral case we get no free parameters: accuracy is maximal.

(operator = DD[4] /. AccuracySolution[4, 4, 1]) // 
      MatrixForm

check[operator, 4, 4, 1] // MatrixForm

If we lower the accuracy keeping the number of points we get free parameters:

(operator = DD[4] /. AccuracySolution[4, 3, 1]) // 
      MatrixForm

check[operator, 4, 3, 1] // MatrixForm

I.4.b. Examples: second-order systems

(operator = DD[4] /. AccuracySolution[4, 2, 2]) // 
      MatrixForm

check[operator, 4, 2, 2] // MatrixForm

(operator = DD[7] /. AccuracySolution[7, 6, 2]) // 
      MatrixForm

check[operator, 7, 6, 2] // MatrixForm

I.4.c. Examples: higher order derivatives

This is a fourth derivative with accuracy 6. We need at least 10 points:

(operator = DD[9] /. AccuracySolution[9, 6, 4]) // 
      MatrixForm

check[operator, 9, 6, 4] // MatrixForm

I.5. Elements far from the diagonal

This function gives the positions of the points on the upper-right corner of the boundary part of the operators. We count from 1, and not from 0.

corner[n_, 1] := {{1, n}}
corner[n_, m_] := 
      Union[corner[n, m - 1], Map[({1, 0} + #) &, corner[n, m - 1]], 
        Map[({0, -1} + #) &, corner[n, m - 1]]]

corner[8, 3]

II.A. SBP conditions on first-order operators

Now we closely follow Strand's paper. (What follows is not self-contained. See that paper.) We define a boundary operator of width r:

B1[r_] := 
      Table[If[i == j == 1, -1/2, 0], {i, r}, {j, 
          r}]

B1[5] // MatrixForm

SetAttributes[b, NHoldAll];
b[i_, i_] := 0
b[i_, j_] := -b[j, i] /; i > j
B2[r_] := Table[b[i - 1, j - 1], {i, r}, {j, r}]

B2[5] // MatrixForm

Power-law vectors of length r or s and exponent j. Boundary part:

e[r_, -1] := Table[0, {r}]
e[r_, j_] := (-1)^j Table[i^j, {i, r, 1, -1}]

e[5, 3] // MatrixForm

and interior part:

f[s_, 0] := Table[1, {s}]
f[s_, j_] := Table[i^j, {i, 0, s - 1}]

f[4, 3] // MatrixForm

Cross term:

CC[r_, s_] := 
      Table[If[j - i + r > s, 0, 1] alpha[j - i + r, s], {i, r}, {j, 
          s}]

CC[5, 2] // MatrixForm

Finally, the inner product:

SetAttributes[h, NHoldAll];
h[i_, j_] := h[j, i] /; i > j
H[r_] := Table[h[i - 1, j - 1], {i, r}, {j, r}]

H[5] // MatrixForm

Condition to restrict to diagonal and restricted-full norms:

diagonalH = h[i_, j_] :> 0 /; i =!= j

restrictedH = h[0, _?Positive] -> 0

H[3] /. diagonalH // MatrixForm

H[3] /. restrictedH // MatrixForm

Following Strand we construct the orthogonal vectors:

g[r_, s_, j_] := 
      j H[r] . e[r, j - 1] - B1[r] . e[r, j] - CC[r, s] . f[s, j]

g[5, 3, 2] // MatrixForm

Compatibility conditions for the H[i,j]. We have (tau+1)(tau+2)/2 conditions (actually one less because the first one is always trivially obeyed) for r(r+1)/2 elements. It is clear that we need r>=tau+1. If we want to have a diagonal norm, we need to impose additional conditions on the H[i,j]. That means that  we need a larger r.

comp[r_, s_, tau_] := 
      Flatten@Table[
          e[r, i] . g[r, s, j] + e[r, j] . g[r, s, i] == 0, {i, 0, 
            tau}, {j, 0, i}]

comp[2, 1, 1] // Simplify // ColumnForm

Useful shortcuts:

bvars[r_] := 
      DeleteCases[Union[Flatten[B2[r]] /. -x_b -> x], 0]

bvars[3]

hvars[r_] := Union@Flatten[H[r]]

hvars[3]

xs[n_] := 
      Table[ToExpression["x" <> ToString[i]], {i, n}]

xs[7]

This function organizes a system of equations:

FlattenEquations[Equal[l1_List, l2_List]] := 
      Flatten[FlattenEquations /@ 
          Thread[Equal[l1, l2]]]
FlattenEquations[x_List] := 
      Flatten[FlattenEquations /@ x]
FlattenEquations[x_] := x

FlattenEquations[{a, b, c} == {1, 2, 3}]

Now we construct the equations to be imposed:

Equations[r_, s_, tau_] := 
      FlattenEquations@
        Table[B2[r] . e[r, j] == g[r, s, j], {j, 0, tau}]

Construct full operator. Arguments:
	r: number of points on the boundary
	s: the accuracy order of the interior stencil is 2s
	tau: accuracy order of the boundary derivative. We need r>=tau+1.

interior[s_] := 
      Table[alpha[nu - i, s], {i, 0, 2  s - 1}, {nu, 0, 2  s - 1}]

zeros[m_, n_] := Table[0, {m}, {n}]

BuildD1[r_, s_, tau_, restrict_: {}] := 
      Module[{matrix = zeros[r + 2  s, r + 2  s], 
          inverseH = 
            Inverse[H[r] /. restrict]}, 
matrix[[Range[
                r], Range[r]]] = 
          inverseH . (B1[r] + B2[r]) /. 
            restrict; 
matrix[[Range[r], 
              Range[r + 1, r + s]]] = 
          inverseH . CC[r, s] /. 
            restrict; 
matrix[[Range[r + 1, r + s], 
              Range[r]]] = -Transpose[CC[r, s]] /. 
            restrict; 
matrix[[Range[r + 1, r + 2  s], 
              Range[r + 1, r + 2  s]]] = 
          interior[
            s]; 
matrix
]

Note that we do not include a Simplify command at the end, because it could take very long. The user must decide whether to use it or not.

BuildD1[2, 1, 1] // MatrixForm

Check only the first r+s points:

checkder1[r_, s_, tau_, other___] := 
      Drop[check[der1[r, s, tau, other], r + 2  s - 1, 2  s, 
          1], -s]

This is the object studied by Peter, directly obtained from the previous function:

ABTE2der1[r_, s_, tau_, other___] := 
Simplify[
        Plus @@ Power[
              Take[(Last /@ checkder1[r, s, tau, other]), r] /. 
                HoldPattern[Derivative[_][f][x]] -> 1, 2]/
          r]

II.B. Results for first-order operators

Given tau>=1 and s>=1 (the order at the boundary is tau and the order in the interior is 2s), we need r>=tau+1 and r>=s. In general there are some free parameters that we can use to get particular forms of H, or improve the previous inequalities. Assume tau+1<=2s (that is, we do not require higher orders at the boundary than in the interior). The best we can get is:
	full norm: r=tau+1, 2s=tau+1 if tau odd, 2s=tau+2 if tau even
	restricted full norm: r=tau+2, 2s=tau+1 if tau odd, 2s=tau+2 if tau even
	diagonal norm:	r=2tau, r>=s>=tau

The equations must be solved in this order:
    	- Solve the compatibility conditions for H, perhaps getting a number of free parameters x1, x2, ... in the nondiagonal cases.
    	- Solve the equations for B, perhaps getting more free parameters.
    	- Impose conditions to fix those free parameters.
    	
Note that full and restricted full norms are slow to compute for r>6.

Full H, r=2, s=1, tau=1. Diagonal in the MB case

Solve for the inner product. There is a free parameter:

parameters = {h[0, 1] -> x1};

Hsol = Solve[comp[2, 1, 1] /. parameters, hvars[2]][[1]]~Join~
        parameters

(norm[2, 1, 1] = H[2] /. Hsol) // MatrixForm

Solve for the derivative operator. There is no free parameter:

Bsol = Solve[Equations[2, 1, 1] /. Hsol, 
          bvars[2]][[1]]

All solutions together:

sol = Join[Hsol, Bsol]

(der1[2, 1, 1] = Simplify@BuildD1[2, 1, 1, sol]) // 
      MatrixForm

For minimum bandwidth choose the corner element to be zero:

solMB = Solve[der1[2, 1, 1][[1, 3]] == 0, x1][[1, 
          1]]

(norm[2, 1, 1, MB] = norm[2, 1, 1] /. solMB) // 
      MatrixForm

(der1[2, 1, 1, MB] = der1[2, 1, 1] /. solMB) // 
      MatrixForm

Checks:

checkder1[2, 1, 1] // MatrixForm

ABTE2 = ABTE2der1[2, 1, 1] /. h^2 -> 1

Plot[Evaluate[ABTE2], {x1, -3, 3}, 
      PlotRange -> {0, 2}]

{junk, rule} = FindMinimum[ABTE2, {x1, 0}]

checkder1[2, 1, 1] /. rule

It is interesting to see that even though the following operator has r=2, it effectively has r=1, what makes the ABTE nonlinear in h:

checkder1[2, 1, 1, MB] // MatrixForm

ABTE2der1[2, 1, 1, MB] // Expand

Sqrt[1/8]

Comparing with Peter we get c1=1/2, which is right, but we do not get the same ABTE.

Full H, r=2, s=2, tau=1. Diagonal in the MB case

Solve for the inner product. There is a free parameter:

parameters = {h[0, 1] -> x1};

Hsol = Solve[comp[2, 2, 1] /. parameters, hvars[2]][[1]]~Join~
        parameters

(norm[2, 2, 1] = H[2] /. Hsol) // MatrixForm

Solve for the derivative operator. There is no free parameter:

Bsol = Solve[Equations[2, 2, 1] /. Hsol, 
          bvars[2]][[1]]

All solutions together:

sol = Join[Hsol, Bsol]

(der1[2, 2, 1] = Simplify@BuildD1[2, 2, 1, sol]) // 
      MatrixForm

For minimum bandwidth choose the corner element to be zero:

solMB = Solve[der1[2, 2, 1][[1, 4]] == 0, x1][[1, 
          1]]

(norm[2, 2, 1, MB] = norm[2, 2, 1] /. solMB) // 
      MatrixForm

(der1[2, 2, 1, MB] = der1[2, 2, 1] /. solMB) // 
      MatrixForm

Checks:

checkder1[2, 2, 1, MB] // MatrixForm

It seems that a better option could be x1=1/8 because we get higher accuracy at the boundary. This is analyzed in the case (2,2,2) next.

checkder1[2, 2, 1] // MatrixForm

ABTE2 = ABTE2der1[2, 2, 1] /. h^2 -> 1

{junk, rule} = FindMinimum[ABTE2, {x1, 0}]

checkder1[2, 2, 1] /. x1 -> 1/8 // MatrixForm

Plot[Evaluate[ABTE2], {x1, -1, 1}, 
      PlotRange -> {0, 2}]

Full H, r=2, s=2, tau=2. No free parameters

There are no free parameters now:

Hsol = Solve[comp[2, 2, 2], hvars[2]][[1]]

(norm[2, 2, 2] = H[2] /. Hsol) // MatrixForm

Solve for the derivative operator. There is no free parameter:

Bsol = Solve[Equations[2, 2, 2] /. Hsol, 
          bvars[2]][[1]]

All solutions together:

sol = Join[Hsol, Bsol]

(der1[2, 2, 2] = BuildD1[2, 2, 2, sol]) // MatrixForm

Checks:

checkder1[2, 2, 2] // MatrixForm

Diagonal H, r=4, s=2, tau=2. No free parameters

Solve for the inner product:

Hsol = Solve[comp[4, 2, 2] /. diagonalH, hvars[4]][[1]]~
        Join~{diagonalH}

(norm[4, 2, 2, diag] = H[4] /. Hsol) // MatrixForm

Solve for the derivative operator. There is no free parameter:

Bsol = Solve[Equations[4, 2, 2] /. Hsol, 
          bvars[6]][[1]]

All solutions together:

sol = Join[Hsol, Bsol]

(der1[4, 2, 2, diag] = BuildD1[4, 2, 2, sol]) // 
      MatrixForm

Checks:

checkder1[4, 2, 2, diag] // MatrixForm

% // N // MatrixForm

ABTE2der1[4, 2, 2, diag] /. h^4 -> 1

% // N

Sqrt[%]

All numbers coincide with Peter's.

Diagonal H, r=6, s=3, tau=3. One free parameter

Solve for the inner product:

Hsol = Solve[comp[6, 3, 3] /. diagonalH, hvars[6]][[1]]~
        Join~{diagonalH}

(norm[6, 3, 3, diag] = H[6] /. Hsol) // MatrixForm

Solve for the derivative operator. There is a single free parameter:

parameters = {b[4, 5] -> x1};

Bsol = Solve[Equations[6, 3, 3] /. diagonalH /. Hsol /. 
              parameters, bvars[6]][[1]]~Join~parameters

All solutions together:

sol = Join[Hsol, Bsol]

(der1[6, 3, 3, diag] = BuildD1[6, 3, 3, sol]) // 
      MatrixForm

For minimum bandwidth choose the corner element to be zero:

solMB = Solve[der1[6, 3, 3, diag][[1, 6]] == 0, 
          x1][[1, 1]]

norm[6, 3, 3, diag, MB] = norm[6, 3, 3, diag];

(der1[6, 3, 3, diag, MB] = der1[6, 3, 3, diag] /. solMB) // 
      MatrixForm

Checks:

checkder1[6, 3, 3, diag, MB] // MatrixForm

% // N // MatrixForm

checkder1[6, 3, 3, diag] // MatrixForm

ABTE2 = ABTE2der1[6, 3, 3, diag] /. h^6 -> 1

{junk, rule} = FindMinimum[ABTE2, {x1, 0}]

checkder1[6, 3, 3, diag] /. rule // MatrixForm

ABTE2 /. rule

Sqrt[%]

All numbers coincide with Peter's.

Diagonal H, r=8, s=4, tau=4. Three free parameters

Solve for the inner product:

Hsol = Solve[comp[8, 4, 4] /. diagonalH, hvars[8]][[1]]~
        Join~{diagonalH}

(norm[8, 4, 4, diag] = H[8] /. Hsol) // MatrixForm

Solve for the derivative operator. There are three free parameters:

parameters = {b[5, 6] -> x1, b[5, 7] -> x2, 
          b[6, 7] -> x3};

Bsol = Solve[Equations[8, 4, 4] /. Hsol /. parameters, 
            bvars[8]][[1]]~Join~parameters

All solutions together:

sol = Join[Hsol, Bsol];

(der1[8, 4, 4, diag] = BuildD1[8, 4, 4, sol]) // 
      MatrixForm

For minimum bandwidth choose the corner element to be zero:

solMB = Solve[{der1[8, 4, 4, diag][[1, 7]] == 
              0, der1[8, 4, 4, diag][[1, 8]] == 
              0, der1[8, 4, 4, diag][[2, 8]] == 0}, {x1, x2, 
            x3}][[1]]

norm[8, 4, 4, diag, MB] = norm[8, 4, 4, diag];

(der1[8, 4, 4, diag, MB] = der1[8, 4, 4, diag] /. solMB) // 
      MatrixForm

Checks:

First, the diagonal and MB case:

checkder1[8, 4, 4, diag, MB] // MatrixForm

% // N // MatrixForm

Second, the diagonal case:

checkder1[8, 4, 4, diag] // MatrixForm

ABTE2 = ABTE2der1[8, 4, 4, diag] /. h^8 -> 1

{junk, rules} = 
      FindMinimum[ABTE2, {{x1, 0}, {x2, 0}, {x3, 0}}, 
        AccuracyGoal -> 16, WorkingPrecision -> 20]

checkder1[8, 4, 4, diag] /. rules // MatrixForm

ABTE2 /. rules

Sqrt[%]

All numbers coincide with Peter's.

I don't remember where these numbers come from:

checkder1[8, 4, 4, diag] /. {x1 -> 0.541`20, 
          x2 -> -0.0675`20, x3 -> 0.748`20} // 
      MatrixForm

Now we try to get the norm given by MN:

MNsol = Solve[{der1[8, 4, 4, diag][[1, 2]] == 
              5544277/5992556, der1[8, 4, 4, diag][[1, 3]] == 
              198794991/
                29962780, der1[8, 4, 4, diag][[3, 2]] == 
              9444709/415220}, {x1, x2, x3}][[1]]

(der1[8, 4, 4, diag, MN] = 
          der1[8, 4, 4, diag] /. MNsol // Simplify) // 
      MatrixForm

norm[8, 4, 4, diag, MN] = norm[8, 4, 4, diag];

Checks:

checkder1[8, 4, 4, diag, MN] // MatrixForm

% // N // MatrixForm

Here we see what MN get with their choice of parameters: there is accuracy h^6 and h^5 at the final points of the boundary, rather than h^4.

Diagonal H, r=10, s=5, tau=5. The norm is not positive definite

Solve for the inner product:

Hsol = 
      Join[Reduce[comp[10, 5, 5] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[10, 5, 5, diag] = H[10] /. Hsol) // MatrixForm

Note that the norm is not definite positive. Therefore this operator cannot be used.

Diagonal H, r=11, s=5, tau=5. No preferred ABTE minimum

Solve for the inner product. Note that tau=6 is not possible:

Reduce[comp[11, 5, 6] /. diagonalH]

Hsol = 
      Join[Reduce[comp[11, 5, 5] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[11, 5, 5, diag] = H[11] /. Hsol /. h[10, 10] -> X) // 
      MatrixForm

The range for X is extremely narrow:

Xcondition = 
      Reduce[Map[# > 0 &, 
          Transpose[norm[11, 5, 5, diag], {1, 1}]]]

% // N

X1 = First[Xcondition];
X2 = Last[Xcondition];

The eigenvalues of the norm are very sensitive to X:

TableForm[Transpose[
              norm[11, 5, 5, diag], {1, 
                1}] /. {{X -> X1}, {X -> X2}} // Transpose // N, 
      TableHeadings -> {None, {"At X1", "At X2"}}]

Solve for the derivative operator. There are 10+1 free parameters:

parameters = {h[10, 10] -> X, b[6, 7] -> x1, 
          b[6, 8] -> x2, b[6, 9] -> x3, b[6, 10] -> x4, 
          b[7, 8] -> x5, b[7, 9] -> x6, b[7, 10] -> x7, 
          b[8, 9] -> x8, b[8, 10] -> x9, 
          b[9, 10] -> x10};

Bsol = Solve[Equations[11, 5, 5] /. Hsol /. parameters, 
            bvars[11]][[1]]~Join~parameters

All solutions together:

sol = Join[Hsol, Bsol];

(der1[11, 5, 5, diag] = Simplify@BuildD1[11, 5, 5, sol]) // 
      MatrixForm

der = der1[11, 5, 5, diag];

For minimum bandwidth choose the corner element to be zero:

solMB = 
      First@Solve[der[[Sequence @@ #]] == 0 & /@ corner[11, 4], 
          xs[10]]

norm[11, 5, 5, diag, MB] = norm[11, 5, 5, diag];

(der1[11, 5, 5, diag, MB] = 
          Simplify[der1[11, 5, 5, diag] /. solMB]) // MatrixForm

Checks:

First, the diagonal and MB case:

checkder1[11, 5, 5, diag, MB] // MatrixForm

Second, the diagonal case:

checkder1[11, 5, 5, diag] // MatrixForm

Now we want to minimize

ABTE2 = ABTE2der1[11, 5, 5, diag] /. h^10 -> 1

None of the denominators vanish in the range for X, but two of them vanish on the boundaries:

denominators = 
      Union[Denominator /@ List @@ ABTE2] /. _.  Power[x_, 2] -> 
          x

Flatten[Union[Solve[# == 0, X] & /@ denominators], 
      1]

Xcondition /. %

denominators /. X->X1

denominators /. X -> X2

For example:

Thread[Rule[xs[10], 0]]

ABTE20 = ABTE2 /. % // Simplify

Plot[Log[10, ABTE20], {X, X1, X2}, Frame -> True]

Peter's choice:

{junk, rules} = 
      FindMinimum[
        ABTE2 /. X -> 101/100, {{x1, 0}, {x2, 0}, {x3, 0}, {x4, 0}, {x5, 
            0}, {x6, 0}, {x7, 0}, {x8, 0}, {x9, 0}, {x10, 0}}, 
        WorkingPrecision -> 20, PrecisionGoal -> 16, 
        AccuracyGoal -> 16]

Sqrt[junk]

checkder1[11, 5, 5, diag] /. X -> 101/100 /. rules // 
      MatrixForm

All numbers coincide with Peter's.

This function finds the minimum of ABTE for a given value of X. The factor k,between 0 and 1, is a linear relative scale between X1 and X2, respectively.

FMABTE2[k_] := 
      FindMinimum[
        ABTE2 /. X -> (1 - k) X1 + k X2, {{x1, 0}, {x2, 0}, {x3, 
            0}, {x4, 0}, {x5, 0}, {x6, 0}, {x7, 0}, {x8, 0}, {x9, 0}, {x10, 
            0}}, WorkingPrecision -> 20, PrecisionGoal -> 16, 
        AccuracyGoal -> 16]

Create a set of points with logarithmic accummulation at 0 and 1:

powers2 = NestList[#/2 &, 1, 10]

ks = 
      Join[Drop[Reverse[powers2], -1], Drop[1 - powers2, 2]]

Sqrt[0.38292897091503204739335557883440543568`20. ]

ABTE2s = First /@ FMABTE2 /@ ks

ListPlot[Transpose[{ks, ABTE2s}], PlotJoined -> True]

It is clear that the minimum is found on the left boundary, but that gives a norm which is not positive definite. The minimum value of ABTE2 does not change much with X anyway: from 0.38 to 0.59. Note that we cannot know that we are actually finding the absolute minimum

Diagonal H, r=11, s=6, tau=5. No preferred ABTE minimum

Hsol = 
      Join[Reduce[comp[11, 6, 5] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[11, 6, 5, diag] = H[11] /. Hsol /. h[10, 10] -> X) // 
      MatrixForm

The range for X is extremely narrow:

Xcondition = 
      Reduce[Map[# > 0 &, 
          Transpose[norm[11, 6, 5, diag], {1, 1}]]]

% // N

X1 = First[Xcondition];
X2 = Last[Xcondition];

Solve for the derivative operator. There are three free parameters:

parameters = {h[10, 10] -> X, b[6, 7] -> x1, 
          b[6, 8] -> x2, b[6, 9] -> x3, b[6, 10] -> x4, 
          b[7, 8] -> x5, b[7, 9] -> x6, b[7, 10] -> x7, 
          b[8, 9] -> x8, b[8, 10] -> x9, 
          b[9, 10] -> x10};

Bsol = Solve[Equations[11, 6, 5] /. Hsol /. parameters, 
            bvars[11]][[1]]~Join~parameters

All solutions together:

sol = Join[Hsol, Bsol];

(der1[11, 6, 5, diag] = Simplify@BuildD1[11, 6, 5, sol]) // 
      MatrixForm

der = der1[11, 6, 5, diag];

For minimum bandwidth choose the corner element to be zero:

solMB = 
      First@Solve[der[[Sequence @@ #]] == 0 & /@ corner[11, 4], 
          xs[10]]

norm[11, 6, 5, diag, MB] = norm[11, 6, 5, diag];

(der1[11, 6, 5, diag, MB] = 
          Simplify[der1[11, 6, 5, diag] /. solMB]) // MatrixForm

Checks:

First, the diagonal and MB case:

checkder1[11, 6, 5, diag, MB] // MatrixForm

Second, the diagonal case:

checkder1[11, 6, 5, diag] // MatrixForm

ABTE2 = ABTE2der1[11, 6, 5, diag] /. h^10 -> 1

None of the denominators vanish in the range for X:

denominators = 
      Union[Denominator /@ List @@ ABTE2[[2]]] /. _.  
            Power[x_, 2] -> x

Flatten[Union[Solve[# == 0, X] & /@ denominators], 
      1]

Xcondition /. %

Some denominators vanish at X1 and X2. That will make Mathematica fail because ABTE will tend to explode:

denominators /. X->X1

denominators /. X -> X2

Thread[Rule[xs[10], 0]]

ABTE20 = ABTE2 /. % // Simplify

Plot[Log[10, ABTE20], {X, X1, X2}, Frame -> True]

This function finds the minimum of ABTE for a given value of X. The factor k,between 0 and 1, is a linear relative scale between X1 and X2, respectively.

FMABTE2[k_] := 
      FindMinimum[
        ABTE2 /. X -> (1 - k) X1 + k X2, {{x1, 0}, {x2, 0}, {x3, 
            0}, {x4, 0}, {x5, 0}, {x6, 0}, {x7, 0}, {x8, 0}, {x9, 0}, {x10, 
            0}}, WorkingPrecision -> 20, PrecisionGoal -> 16, 
        AccuracyGoal -> 16]

Create a set of points with logarithmic accummulation at 0 and 1:

powers2 = NestList[#/2 &, 1, 10]

ks = 
      Join[Drop[Reverse[powers2], -1], Drop[1 - powers2, 2]]

ABTE2s = First /@ FMABTE2 /@ ks

ListPlot[Transpose[{ks, ABTE2s}], PlotJoined -> True]

Diagonal H, r=12, s=6, tau=5. Nothing useful

Hsol = 
      Join[Reduce[comp[12, 6, 5] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[12, 6, 5, 
            diag] = H[12] /. Hsol /. {h[10, 10] -> X, 
              h[11, 11] -> Y}) // MatrixForm

We find just a very narrow range of accepted values:

Xconditions = 
      Reduce[Map[# > 0 &, 
          Transpose[norm[12, 6, 5, diag], {1, 1}]]]

plot[
        HoldPattern[
          And[Inequality[YL_, _, Y, _, YR_], 
            Inequality[fDown_, _, X, _, fUp_]]]] := 
      Plot[{fDown, fUp}, {Y, YL, YR}, 
        DisplayFunction -> Identity]

Show[plot /@ List @@ Xconditions, 
      DisplayFunction -> $DisplayFunction, Frame -> True, 
      Axes -> False]

Diagonal H, r=12, s=6, tau=6. The norm is not positive definite

Hsol = 
      Join[Reduce[comp[12, 6, 6] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[12, 6, 6, diag] = H[12] /. Hsol) // MatrixForm

Diagonal H, r=13, s=6, tau=6. The norm is not positive definite

Hsol = 
      Join[Reduce[comp[13, 6, 6] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[13, 6, 6, diag] = H[13] /. Hsol /. h[12, 12] -> X) // 
      MatrixForm

There is no positive definite case:

Xcondition = 
      Reduce[Map[# > 0 &, 
          Transpose[norm[13, 6, 6, diag], {1, 1}]]]

Diagonal H, r=14, s=7, tau=7. The norm is not positive definite

Hsol = 
      Join[Reduce[comp[14, 7, 7] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[14, 7, 7, diag] = H[14] /. Hsol) // MatrixForm

Diagonal H, r=15, s=7, tau=7. The norm is not positive definite

Hsol = 
      Join[Reduce[comp[15, 7, 7] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[15, 7, 7, diag] = H[15] /. Hsol /. h[14, 14] -> X) // 
      MatrixForm

There is no positive definite case:

Xcondition = 
      Reduce[Map[# > 0 &, 
          Transpose[norm[15, 7, 7, diag], {1, 1}]]]

Diagonal H, r=16, s=8, tau=8. The norm is not positive definite

Hsol = 
      Join[Reduce[comp[16, 8, 8] /. diagonalH] // 
          ToRules, {diagonalH}]

(norm[16, 8, 8, diag] = H[16] /. Hsol) // MatrixForm

Full H, r=3, s=2, tau=2. Two free parameters

Solve for the inner product. There are two free parameters:

parameters = {h[0, 2] -> x1, 
          h[1, 2] -> x2 - 1/16};

Hsol = Solve[comp[3, 2, 2] /. parameters, hvars[3]][[1]]~Join~
        parameters

(HH = norm[3, 2, 2] = H[3] /. Hsol) // MatrixForm

Solve for the derivative operator. There are no additional free parameters:

Bsol = Solve[Equations[3, 2, 2] /. Hsol, 
          bvars[3]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[3, 2, 2] = Simplify@BuildD1[3, 2, 2, sol]) // 
      MatrixForm

ff = Denominator[DD[[1, 1]]]

ff DD[[{1, 2, 3}]] // Simplify // Expand // 
      MatrixForm

We can use x1 and x2 to minimize the bandwidth. There are three solutions:

sols = 
      Solve[{%[[1, 5]] == 0, %[[2, 5]] == 0}, {x1, 
          x2}]

Only the first one gives a positive definite HH:

(norm[3, 2, 2, MB] = HH /. sols[[1]] // Simplify) // 
      MatrixForm

Eigenvalues[%]

% // N

Eigenvalues[HH /. sols[[2]]] // N

Eigenvalues[HH /. sols[[3]]] // N

Finally:

(der1[3, 2, 2, MB] = DD /. sols[[1]]) // MatrixForm

Checks:

checkder1[3, 2, 2] // MatrixForm

checkder1[3, 2, 2, MB] // MatrixForm

ABTE2der1[3, 2, 2, MB]

Sqrt[%] /. h^8 -> 0 /. h^4 -> 1 // N

{junk, rules} = 
      FindMinimum[
        ABTE2der1[3, 2, 2] /. h -> 1, {{x1, 0}, {x2, 0}, {x3, 0}}, 
        WorkingPrecision -> 20, PrecisionGoal -> 16, 
        AccuracyGoal -> 16]

checkder1[3, 2, 2] /. rules // MatrixForm

Full H, r=3, s=3, tau=2. Two free parameters

Solve for the inner product. There are two free parameters:

parameters = {h[0, 2] -> x1, h[1, 2] -> x2};

Hsol = Solve[comp[3, 3, 2] /. parameters, hvars[3]][[1]]~Join~
        parameters

(HH = norm[3, 3, 2] = H[3] /. Hsol) // MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[3, 3, 2] /. Hsol, 
          bvars[3]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[3, 3, 2] = Simplify@BuildD1[3, 3, 2, sol]) // 
      MatrixForm

ff = Denominator[DD[[1, 1]]]

ff DD[[{1, 2, 3}]] // Simplify // Expand // 
      MatrixForm

We can use x1 and x2 to minimize the bandwidth. There are three solutions:

sols = 
      Solve[{%[[1, 6]] == 0, %[[2, 6]] == 0}, {x1, 
          x2}]

Only the first one gives a positive definite HH:

(norm[3, 3, 2, MB] = HH /. sols[[1]] // Simplify) // 
      MatrixForm

Eigenvalues[%]

% // N

Eigenvalues[HH /. sols[[2]]] // N

Eigenvalues[HH /. sols[[3]]] // N

Finally:

(der1[3, 3, 2, MB] = DD /. sols[[1]]) // MatrixForm

Checks:

checkder1[3, 3, 2] // MatrixForm

checkder1[3, 3, 2, MB] // MatrixForm

ABTE2der1[3, 3, 2, MB]

Sqrt[%] /. h^4 -> 1 // N

{junk, rules} = 
      FindMinimum[ABTE2der1[3, 3, 2] /. h^4 -> 1, {{x1, 0}, {x2, 0}}, 
        WorkingPrecision -> 20, PrecisionGoal -> 16, 
        AccuracyGoal -> 16]

checkder1[3, 3, 2] /. rules // MatrixForm

Full H, r=4, s=2, tau=3. Two free parameters

Solve for the inner product. There are two free parameters:

parameters = {h[0, 0] -> -4421/54 + x1/12 + x2/36, 
          h[0, 1] -> 235099/864 - 3  x1/8 - x2/12};

Hsol = Solve[comp[4, 2, 3] /. parameters, hvars[4]][[1]]~Join~
        parameters

(HH = norm[4, 2, 3] = H[4] /. Hsol) // MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[4, 2, 3] /. Hsol, 
          bvars[4]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[4, 2, 3] = Simplify@BuildD1[4, 2, 3, sol]) // 
      MatrixForm

ff = Denominator[DD[[1, 1]]]

ff DD[[{1, 2, 3, 4}]] // Simplify // Expand // 
      MatrixForm

We can use x1 and x2 to minimize the bandwidth. There are four solutions:

sols = 
      Solve[{%[[1, 6]] == 0, %[[2, 6]] == 0}, {x1, 
          x2}]

Only the second one gives a positive definite HH:

(norm[4, 2, 3, MB] = HH /. sols[[2]] // Simplify) // 
      MatrixForm

Eigenvalues[%] // N

Eigenvalues[HH /. sols[[1]]] // N

Eigenvalues[HH /. sols[[3]]] // N

Eigenvalues[HH /. sols[[4]]] // N

Finally:

(der1[4, 2, 3, MB] = DD /. sols[[2]] // Simplify) // 
      MatrixForm

Checks:

checkder1[4, 2, 3] // MatrixForm

checkder1[4, 2, 3, MB] // MatrixForm

% // N // MatrixForm

ABTE2der1[4, 2, 3, MB] /. h^6 -> 1 // N

Note that the ABTE has this form. It is almost constant with respect to x2, but decreases towards -Infinity:

Plot3D[
      Evaluate[ABTE2der1[4, 2, 3] /. h -> 1], {x1, 180, 
        220}, {x2, -1000000, 0}]

FMABTE2[y_] := 
      FindMinimum[ABTE2der1[4, 2, 3] /. h^6 -> 1 /. 
          x2 -> y, {x1, 0}, WorkingPrecision -> 30, 
        PrecisionGoal -> 20, AccuracyGoal -> 20, 
        MaxIterations -> 1000]

FMABTE2 /@ Table[-10^i, {i, 0, 10}]

Now we try to get the norm given by MN:

r1 = -(2177  Sqrt[295369] - 1166427)/
          25488;
r2 = (66195  Sqrt[53] Sqrt[5573] - 35909375)/
          101952;

MNsol = Solve[{norm[4, 2, 3][[1, 1]], norm[4, 2, 3][[1, 
                  2]]} == {-(216  r2 + 2160  r1 - 2125)/
                12960, (81  r2 + 675  r1 + 415)/540}, {x1, 
            x2}][[1]]

(norm[4, 2, 3, MN] = norm[4, 2, 3] /. MNsol // Simplify) // 
      MatrixForm

Eigenvalues[%] // N

We can now check other values:

-(72  r2 + 720  r1 + 445)/1440 - norm[4, 2, 3][[1, 3]] /. 
        MNsol // Simplify

which gives

(der1[4, 2, 3, MN] = der1[4, 2, 3] /. MNsol // Simplify) // N // 
      MatrixForm

instead of

der1[4, 2, 3, MB] // Simplify // N // MatrixForm

Full H, r=4, s=4, tau=3

Solve for the inner product. There are two free parameters:

parameters = {h[0, 3] -> x1, h[1, 3] -> x2};

Hsol = Solve[comp[4, 4, 3] /. parameters, hvars[4]][[1]]~Join~
        parameters

(HH = norm[4, 4, 3] = H[4] /. Hsol) // MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[4, 4, 3] /. Hsol, 
          bvars[4]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[4, 4, 3] = 
            Simplify@BuildD1[4, 4, 3, sol] // Simplify) // 
      MatrixForm

ff = Denominator[DD[[1, 1]]]

ff DD[[{1, 2, 3, 4}]] // Simplify // Expand // 
      MatrixForm

We can use x1 and x2 to minimize the bandwidth. There are four solutions:

sols = 
      Solve[{%[[1, 8]] == 0, %[[2, 8]] == 0}, {x1, x2}] // 
        Simplify

Only the fourth one gives a positive definite HH:

(norm[4, 4, 3, MB] = HH /. sols[[4]] // Simplify) // 
      MatrixForm

Eigenvalues[%] // N

Eigenvalues[HH /. sols[[1]]] // N

Eigenvalues[HH /. sols[[2]]] // N

Eigenvalues[HH /. sols[[3]]] // N

Finally:

(norm[4, 4, 3, MB] = HH /. sols[[4]] // Simplify) // N // 
      MatrixForm

(der1[4, 4, 3, MB] = DD /. sols[[4]] // Simplify) // N // 
      MatrixForm

Checks:

checkder1[4, 4, 3] // MatrixForm

checkder1[4, 4, 3, MB] // MatrixForm

% // N // MatrixForm

Full H, r=6, s=3, tau=5

Solve for the inner product. There are three free parameters:

parameters = {h[5, 5] -> x1, h[4, 5] -> x2, 
          h[3, 5] -> x3};

Hsol = Solve[comp[6, 3, 5] /. parameters, hvars[6]][[1]]~Join~
        parameters

(HH = norm[6, 3, 5] = H[6] /. Hsol) // MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[6, 3, 5] /. Hsol, 
          bvars[6]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[6, 3, 5] = Simplify@BuildD1[6, 3, 5, sol]) // 
      MatrixForm

We can use x1, x2 and x3 to minimize the bandwidth. There are 8 solutions, 2 of which are complex

sols = 
      NSolve[{DD[[1, 9]] == 0, DD[[1, 8]] == 0, 
          DD[[2, 9]] == 0}, {x1, x2, x3}, 
        WorkingPrecision -> 20]

Only solution 6 gives a positive definite HH:

HH /. sols[[1]] // Eigenvalues

HH /. sols[[4]] // Eigenvalues

HH /. sols[[5]] // Eigenvalues

HH /. sols[[6]] // Eigenvalues

HH /. sols[[7]] // Eigenvalues

HH /. sols[[8]] // Eigenvalues

(norm[6, 3, 5, MB] = HH /. sols[[6]]) // MatrixForm

(der1[6, 3, 5, MB] = DD /. sols[[6]]) // MatrixForm

Check:

checkder1[6, 3, 5, MB] // MatrixForm

Full restricted H, r=5, s=2, tau=3. This is Peter's case.

Solve for the inner product. There are three free parameters:

parameters = {h[0, 0] -> x1, h[1, 1] -> x2, 
          h[1, 2] -> x3};

Hsol = Solve[comp[5, 2, 3] /. restrictedH /. parameters, 
            hvars[5]][[1]]~Join~parameters~
        Join~{restrictedH}

(HH = norm[5, 2, 3, restr] = H[5] /. Hsol) // 
      MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[5, 2, 3] /. Hsol, 
          bvars[5]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

(DD = der1[5, 2, 3, restr] = Simplify@BuildD1[5, 2, 3, sol]) // 
      MatrixForm

DD[[1]]

ff = Expand@Denominator[DD[[2, 1]]]

ff DD[[{2, 3, 4, 5}]] // Simplify // Expand // 
      MatrixForm

(ff DD // Simplify // Expand) /. _Plus -> X // 
      MatrixForm

We can minimize the bandwidth with:

sols = 
      Solve[{DD[[1, 5]] == 0, %%[[1, 7]] == 
              0, %%[[2, 7]] == 0}, {x1, x2, x3}] // 
        Simplify

We need to use solution 2:

Eigenvalues[HH /. sols[[1]]] // N

Eigenvalues[HH /. sols[[3]]] // N

Eigenvalues[HH /. sols[[4]]] // N

Eigenvalues[HH /. sols[[2]]] // N

norm[5, 2, 3, restr, 
          MB] = (HH /. sols[[2]] // Simplify);

%[[1]]

%% // N // MatrixForm

der1[5, 2, 3, restr, MB] = 
        DD /. sols[[2]] // Simplify;

%[[1]]

%% // N // MatrixForm

checkder1[5, 2, 3, restr] // MatrixForm

checkder1[5, 2, 3, restr, MB] // N // MatrixForm

Sqrt[ABTE2der1[5, 2, 3, restr, MB]] /. h^6 -> 1 // N

We try to get a smaller value:

{junk90, rules90} = 
      FindMinimum[
        ABTE2der1[5, 2, 3, restr] /. 
          h^6 -> 1, {{x1, 1}, {x2, 2}, {x3, -1}}, 
        WorkingPrecision -> 90, PrecisionGoal -> 50, 
        AccuracyGoal -> 50]

{junk100, rules100} = 
      FindMinimum[
        ABTE2der1[5, 2, 3, restr] /. 
          h^6 -> 1, {{x1, 1}, {x2, 2}, {x3, -1}}, 
        WorkingPrecision -> 100, PrecisionGoal -> 50, 
        AccuracyGoal -> 50]

junk90 - junk100

Read Peter's operator:

<< "~/projects/peter/Restricted_Full_4_3.m";

rule1 = Solve[DD[[1, 1]] == DDEC[[1, 1]], x1][[1, 
          1]]

A check of accuracy:

(DD[[1]] /. rule1) - DDEC[[1]]

(HH[[1]] /. rule1) - HHEC[[1]]

Reduce[{DD[[2, 2]] == DDEC[[2, 2]], 
          DD[[2, 3]] == DDEC[[2, 3]]} /. rule1, {x2, 
        x3}]

rules23 = Solve[%, {x2, x3}]

Max@Abs[(DD /. rules23[[1]] /. rule1) - DDEC]

Max@Abs[(DD /. rules23[[2]] /. rule1) - DDEC]

We need to use rule number 1:

rulesPeter1 = Sort@Join[{rule1}, rules23[[1]]]

ABTE2Peter1 = ABTE2der1[5, 2, 3, restr] /. rulesPeter1 /. 
        h -> 1

The parameters, directly from Peter:

rulesPeter2 = << "~/projects/peter/parameters-4-3.txt"

ABTE2Peter2 = ABTE2der1[5, 2, 3, restr] /. h^6 -> 1 /. 
        rulesPeter2

There are no differences among the numbers supplied by Peter:

(Last /@ rulesPeter1) - (Last /@ rulesPeter2)

ABTE2Peter1 - ABTE2Peter2

The differences in the parameters are of order 10^(-6) and the difference in ABTEs is of order 10^(-11):

(Last /@ rulesPeter2) - (Last /@ rules100)

ABTE2Peter2 - junk90

Sqrt[ABTE2Peter2] - Sqrt[junk90]

Finally we construct the first order operators:

MatrixForm[norm[5, 2, 3, restr, Peter] = HH /. rulesPeter2]

MatrixForm[der1[5, 2, 3, restr, Peter] = DD /. rulesPeter2]

Full restricted H, r=7, s=3, tau=5. Suggested by Peter Diener

Solve for the inner product. There are four free parameters:

parameters = {h[0, 0] -> x1, h[6, 6] -> x2, 
          h[5, 6] -> x3, h[4, 6] -> x4};

Hsol = Solve[comp[7, 3, 5] /. restrictedH /. parameters, 
            hvars[7]][[1]]~Join~parameters~
        Join~{restrictedH}

(HH = norm[7, 3, 5, restr] = H[7] /. Hsol) // 
      MatrixForm

Solve for the derivative operator. There are no additional free parameters:

Bsol = Solve[Equations[7, 3, 5] /. Hsol, 
          bvars[7]][[1]]

All solutions together:

sol = Join[Hsol, Bsol];

This takes very long. Do not use Simplify.

Timing[DD = der1[7, 3, 5, restr] = 
            BuildD1[7, 3, 5, sol];]

DD[[1]]

Example of term in the derivative:

Simplify[DD[[2, 1]]]

We can minimize the bandwidth with:

sols = 
      NSolve[Simplify[{DD[[1, 7]] == 0, DD[[2, 10]] == 0, 
            DD[[2, 9]] == 0, DD[[3, 10]] == 0}], {x1, x2, 
          x3, x4}, WorkingPrecision -> 20]

We need the solution 2:

Eigenvalues[HH /. sols[[1]]] // N

Eigenvalues[HH /. sols[[2]]] // N

Eigenvalues[HH /. sols[[7]]] // N

Eigenvalues[HH /. sols[[8]]] // N

norm[7, 3, 5, restr, 
          MB] = (norm[7, 3, 5, restr] /. sols[[2]] // 
            Simplify);

% // N // MatrixForm

(der1[7, 3, 5, restr, MB] = 
          Rationalize[der1[7, 3, 5, restr] /. sols[[2]]] // Simplify) // 
      MatrixForm

checkder1[7, 3, 5, restr, MB] // MatrixForm

Full restricted H, r=8, s=4, tau=6. Suggested by Peter Diener. This takes too long

Solve for the inner product. There are five free parameters:

parameters = {h[0, 0] -> x1, h[7, 7] -> x2, 
          h[6, 7] -> x3, h[5, 7] -> x4, 
          h[4, 7] -> x5};

Hsol = 
      Join[Solve[comp[8, 4, 6] /. restrictedH, 
              RotateLeft[hvars[8], 1]][[1]] /. parameters, 
        parameters, {restrictedH}]

(HH = norm[8, 4, 6, restr] = H[8] /. Hsol) // 
      MatrixForm

Solve for the derivative operator. There are no free parameters:

Bsol = Solve[Equations[8, 4, 6] /. Hsol, 
          bvars[8]][[1]]

All solutions together:

sol = Dispatch@Join[Hsol, Bsol];

Derivative operator:

First@
      Timing[DD = der1[8, 4, 6, restr] = 
            BuildD1[8, 4, 6, sol]]

First row:

Simplify@DD[[1]]

General structure:

DD /. {_Plus -> X, _Times -> X} // MatrixForm

We can minimize the bandwidth with:

First@
      Timing[eqs = 
          Simplify[{DD[[2, 12]] == 0, DD[[2, 11]] == 0, 
              DD[[3, 11]] == 0, DD[[3, 12]] == 0, 
              DD[[2, 10]] == 0}]]

For some reason Mathematica 5.2 cannot solve this system of equations with 512Mbytes of RAM:

Timing[sols = NSolve[eqs, {x1, x2, x3, x4, x5}]]

Eigenvalues[HH /. sols[[1]]] // N

Eigenvalues[HH /. sols[[2]]] // N

Eigenvalues[HH /. sols[[7]]] // N

Eigenvalues[HH /. sols[[8]]] // N

norm[8, 4, 6, restr, 
          MB] = (norm[8, 4, 6, restr] /. sols[[2]] // 
            Simplify);

% // N // MatrixForm

(der1[8, 4, 6, restr, MB] = 
          Rationalize[der1[8, 4, 6, restr] /. sols[[2]]] // Simplify) // 
      MatrixForm

checkder1[8, 4, 6, restr, MB] // MatrixForm

Dump

TimeUsed[]

MemoryInUse[]

MaxMemoryUsed[]

SetDirectory["~/projects/peter/"]

Exit[]

DumpSave["SBP_LSU.mx"]

Get["SBP_LSU.mx"]

MemoryInUse[]

III.A. SBP conditions on second-order operators

Change of notation from Strand to MN

In order to meet Mattson & Nordstrom (MN) notation we construct "double operators". The translation between Strand to MN operators is done using this function:

DoubleReverse[{}] = {};
DoubleReverse[s_, matrix : {__List}] := 
      Reverse /@ Reverse[s matrix]

der1 to doubleder1:

doubleder1[][r_, s_, tau_, 
          other___] := doubleder1[2  r + 2  s][r, s, tau, 
          other];
doubleder1[total_][r_, s_, tau_, other___] := 
      Module[{matrix = zeros[total, total], der = der1[r, s, tau, other], 
          interior, i}, 
If[total < 2  r, 
          Print["total size is too small."]; 
          Return[]]; 
interior = 
          PadRight[der[[r + 1]], 
            total]; 
matrix[[Range[1, r], 
              Range[1, r + 2  s]]] = 
          der[[Range[r]]]; 
matrix[[Range[
                total - r + 1, total], 
              Range[total - r - 2  s + 1, total]]] = 
          DoubleReverse[-1, der[[Range[r]]]]; 
Do[
          matrix[[i]] = RotateRight[interior, i - (r + 1)], {i, 
            r + 1, total - 
              r}]; 
matrix
]

DOs in Strand notation has been given size r+2s:

der1[3, 2, 2, MB] // MatrixForm

In MN notation they will have size 2r+2s:

doubleder1[][3, 2, 2, MB] // MatrixForm

norm to doublenorm:

doublenorm[][r_, s_, tau_, 
          other___] := doublenorm[2  r + 2  s][r, s, tau, 
          other];
doublenorm[total_][r_, s_, tau_, other___] := 
      Module[{matrix = IdentityMatrix[total], 
          nor = norm[r, s, tau, other]}, 
If[total < 2  r, 
          Print["total size is too small."]; 
          Return[]]; 
matrix[[Range[1, r], 
              Range[1, r]]] = 
          nor; 
matrix[[Range[total - r + 1, total], 
              Range[total - r + 1, total]]] = 
          DoubleReverse[1, 
            nor]; 
matrix
]

Norms in Strand notations are given with size r:

norm[3, 2, 2, MB] // MatrixForm

and now they also have 2r+2s size:

doublenorm[][3, 2, 2, MB] // MatrixForm

Boundary operator:

doubleB[][r_, s_, tau_, other___] := doubleB[2  r + 2  s][r, 
          s, tau, other];
doubleB[total_][r_, s_, tau_, other___] := 
        If[total < 2  r, Print["total size is too small."]; Return[], 
          Table[Which[i == j == 1, -1, 
              i == j == total, 1, True, 0], {i, total}, {j, 
              total}]];

doubleB[][3, 2, 2, MB] // MatrixForm

Check SBP condition for a 1st order operator: HD1 + (HD1)^t = B.Just an example:

With[{Q = doublenorm[][3, 2, 2, MB] . doubleder1[][3, 2, 2, 
              MB]}, Q + Transpose[Q] == doubleB[][3, 2, 2, 
          MB]
]

Construction of second-order operators

First note that a second-order SBP operator is always associated to a particular first-order SBP operator. Therefore we need to supply both parameters describing the first-order and the second-order operators. There is also a second first-order operator S.

Now we construct the matrix A of MN. We assume that the interior stencil is centered and give its half-width s2. That means that there are 2s2+1 points in the stencil and that it is of order 2s2 for the second derivative.

doubleA[][{r_, s_, tau_, other___}, 
          s2_] := doubleA[2  r + 2  s][{r, s, tau, other}, 
          s2];
doubleA[total_][{r_, s_, tau_, other___}, s2_] := 
      Module[{matrix = 
            zeros[total, total], 
interiorstencil = 
            Simplify[-h^2 UFDWeights[2, 2  s2, 
                  s2]]}, 
If[total < 2  r, 
          Print["total size is too small."]; 
          Return[]]; 
matrix[[Range[r], 
              Range[r + s]]] = 
          Table[a[i, j], {i, r}, {j, r + s}]; 
For[i = 1, 
          i ≤ total - 2  r, i++, 
          matrix[[r + i, Range[i + r - s2, i + r + s2]]] = 
            interiorstencil]; 
matrix[[Range[
                total - r + 1, total]]] = 
          DoubleReverse[1, 
            matrix[[Range[r]]]]; 
matrix]

doubleA[][{3, 2, 2}, 1] // MatrixForm

First-derivative on the boundary:

doubleS[size_, rS_] := 
      Module[{matrix = zeros[size, size], 
          list = Table[S[j], {j, rS}]}, 
matrix[[1, 
              Range[1, rS]]] = 
          list; 
matrix[[size, 
              Range[size - rS + 1, size]]] = -Reverse[
              list]; 
matrix]

Remember that we need: rS≥tauS+1

Here we encode the SBP conditions for the a second-order DO. Parameters:
	- First-order parameters:
		- r : number of points on the boundary
		- s1 : semi-width of the first-order DO in the interior
		- tau1 : accuracy degree of the first-order DO on the boundary
		- other : optional flags
	- Second-order parameters:
		- s2 : semi-width of the second-order DO in the interior
		- tau2 : accuracy degree of the second-order DO on the boundary
		- rS : number of points of the boundary stencil of S
		- tauS : accuracy of the first-order DO S

D2conditions[][{r_, s1_, tau1_, other___}, 
          d2__] := D2conditions[2  r + 2  s1][{r, s1, tau1, other}, 
          d2];
D2conditions[total_][{r_, s1_, tau1_, other___}, s2_, tau2_, rS_, 
        tauS_] := Module[{D2, boundary = Range[1, r], 
          interior = Range[r + 1, total - r]}, 
        D2 = Simplify@RootReduce@
                Inverse[doublenorm[total][r, s1, tau1, other], 
                  Method -> 
                    OneStepRowReduction] . (-doubleA[total][{r, s1, 
                      tau1, other}, s2] + doubleB[total][r, s1, tau1] . 
                  doubleS[total, 
                    rS]); 
FlattenEquations@{(doubleS[
                      total, rS][[1]]) . 
                LeftMatrix[total - 1, tauS + 1] == RightMatrix[
                  total - 1, tauS, 1][[1]], 
            D2[[boundary]] . 
                LeftMatrix[total - 1, tau2 + 2] == RightMatrix[total - 1, 
                  tau2, 2][[boundary]]}]

symmetricA := a[j_, i_] :> a[i, j] /; j > i

BuildD2[][{r_, s_, tau_, other___}, 
          d2__] := BuildD2[2  r + 2  s][{r, s, tau, other}, 
          d2];
BuildD2[total_][{r_, s_, tau_, other___}, s2_, tau2_, rS_, tauS_, 
        restrict_: {}] := 
      Simplify@RootReduce[
          Inverse[doublenorm[total][r, s, tau, other], 
                Method -> 
                  OneStepRowReduction] . (-doubleA[total][{r, s, tau, 
                        other}, s2] + doubleB[total][r, s, tau, 
                      other] . doubleS[total, rS]) //. 
            restrict]

Checks:

checkder2[][r_, s_, tau_, other___] := checkder2[2  r + 2  s][
          r, s, tau, other];
checkder2[total_][r_, s_, tau_, other___] := 
      check[der2[total][r, s, tau, other], total - 1, Max[2  s, tau], 
        2]

checkS[operatorS_, size_] := 
      Take[Normal[
          Simplify@
            Series[1/h operatorS[[1]] . 
                  Table[f[x + i h], {i, 0, size - 1}], {h, 0, size}]], 
        2]

Strange bug in Mathematica 5 (5.0, 5.1 and 5.2), already reported to Wolfram.com

Suppose we have

(matrix = doublenorm[4, 2, 3, MB]) // MatrixForm

The matrix is perfectly invertible:

Det[matrix] // Simplify

However Mathematica cannot do it:

Inverse[matrix]

After reporting the bug, Wolfram_Support recommends me changing the standard inversion method to:

inverse = 
      Inverse[matrix,  Method -> OneStepRowReduction] // 
        Simplify

Note the different representations of the same numbers:

(inverse2 = matrix[[{1, 2, 3, 4}, {1, 2, 3, 4}]] // Inverse // 
            Simplify) // MatrixForm

inverse[[{1, 2, 3, 4}, {1, 2, 3, 4}]] // Simplify // 
      MatrixForm

inverse[[{-1, -2, -3, -4}, {-1, -2, -3, -4}]] // Simplify // MatrixForm

This particular problem can be solved by using RootReduce:

(inverse2 = matrix[[{1, 2, 3, 4}, {1, 2, 3, 4}]] // 
                Inverse // RootReduce // Simplify) // 
      MatrixForm

inverse[[{1, 2, 3, 4}, {1, 2, 3, 4}]] // RootReduce // 
        Simplify // MatrixForm

Unfortunately, the result is bad again in Mathematica 5.2 (it worked in 5.0):

inverse[[{-1, -2, -3, -4}, {-1, -2, -3, -4}]] // RootReduce // Simplify // MatrixForm

The problem can be traced back to insufficient internal numerical precision in certain intermediate estimates:

FlattenEquations[
        inverse[[{1, 2, 3, 4}, {1, 2, 3, 4}]] == inverse2] // 
      FullSimplify

eq = (-47759432851319 + 
              114299464067 (Sqrt[174601])) (-28453134961 + 
              68178811 (Sqrt[174601])) == 
        32 (-124358779 + 300529 (Sqrt[174601])) (-340029130647717 + 
              813755314865 (Sqrt[174601]))

We have

N[eq, 16]

FullSimplify[eq]

but

N[eq]

FullForm[eqL = N[eq[[1]]]]

{Precision[eqL], Accuracy[eqL]}

FullForm[eqR = N[eq[[2]]]]

{Precision[eqR], Accuracy[eqR]}

eqL - eqR

%/eqL

This means that the accuracy is five orders of magnitude below what Mathematica reports! Horrible.

III.B. Results for second order-operators

Diagonal H, r=2, s=1, tau=1, s2=1, tau2=1, rS=3, tauS=2. Appendix C.1 of MN

First order operators:

doublenorm[][2, 1, 1, MB] // MatrixForm

doubleder1[][2, 1, 1, MB] // MatrixForm

These are the variables we must find:

allvars = {a[1, 1], a[1, 2], a[1, 3], a[2, 2], a[2, 3], S[1], S[2], 
          S[3]};

and these are the equations they obey:

eqs = D2conditions[][{2, 1, 1, MB}, 1, 1, 3, 2] /. 
        symmetricA

sol = Solve[%, allvars][[1]]~Join~{symmetricA}

(AA = doubleA[][{2, 1, 1, MB}, 1] //. sol) // 
      MatrixForm

There is a zero eigenvalue:

Eigenvalues[AA]

(SS = doubleS[6, 3] /. sol) // MatrixForm

This operator corrects that of MN:

doubleB[][2, 1, 1, MB] . SS // MatrixForm

(DD = der2[6][2, 1, 1, MB] = BuildD2[][{2, 1, 1, MB}, 1, 1, 
              3, 2, sol]) // MatrixForm

checkder2[][2, 1, 1, MB] // MatrixForm

checkS[SS, 6]

Diagonal H, r=4, s=2, tau=2, s2=2, tau2=2, rS=4, tauS=3. Appendix C.2 of MN

First-order operators:

norm[4, 2, 2, diag] // MatrixForm

der1[4, 2, 2, diag] // MatrixForm

Free parameters:

parameters = {a[1, 5] -> y1, a[1, 6] -> y2, 
          a[2, 6] -> y3};

Variables we must solve for:

allvars = 
      Complement[
        Cases[Union[
              Flatten@doubleA[][{4, 2, 2, diag}, 2] /. symmetricA], _a]~
          Join~DeleteCases[
            Union[Flatten@doubleS[12, 4] /. -x_ -> x], 0], 
        First /@ parameters]

Equations they obey:

eqs = D2conditions[][{4, 2, 2, diag}, 2, 2, 4, 3] /. 
            symmetricA /. parameters;

Solution:

sol = Solve[eqs, allvars][[1]]~Join~parameters~
        Join~{symmetricA}

(AA = doubleA[][{4, 2, 2, diag}, 2] //. sol) // 
      MatrixForm

There is always a zero eigenvalue:

Eigenvalues[AA]

Eigenvalues[AA /. {y1 -> 0, y2 -> 0, y3 -> 0}] // 
      N

(SS = doubleS[12, 4] /. sol) // MatrixForm

doubleB[][4, 2, 2, diag] . SS // MatrixForm

Second derivative operator:

(DD = der2[12][4, 2, 2, 
              diag] = BuildD2[][{4, 2, 2, diag}, 2, 2, 4, 3, 
                sol] /. {y1 -> 0, y2 -> 0, y3 -> 0}) // 
      MatrixForm

checkder2[][4, 2, 2, diag] // MatrixForm

checkS[SS, 12]

Diagonal H, r=6, s=3, tau=3, s2=3, tau3=2, rS=5, tauS=4. Appendix C.3 of MN

First-order operators:

norm[6, 3, 3, diag] // MatrixForm

der1[6, 3, 3, diag] // MatrixForm

In order to get MN result we take (don't know the reason)

x1rule = Solve[der1[6, 3, 3, diag][[1, 2]] == 
            104009/54596, x1][[1, 1]]

der1[6, 3, 3, diag] /. x1rule // MatrixForm

Free parameters:

parameters = {a[1, 6] -> y1, a[1, 7] -> y2, 
          a[1, 8] -> y3, a[1, 9] -> y4, a[2, 7] -> y5, 
          a[2, 8] -> y6, a[2, 9] -> y7, a[3, 8] -> y8, 
          a[3, 9] -> y9, a[4, 9] -> y10};

Variables we must solve for:

allvars = 
      Complement[
        Cases[Union[
              Flatten@doubleA[][{6, 3, 3, diag}, 3] /. symmetricA], _a]~
          Join~DeleteCases[
            Union[Flatten@doubleS[18, 5] /. -x_ -> x], 0], 
        First /@ parameters]

Equations they obey:

eqs = 
        Simplify[D2conditions[][{6, 3, 3, diag}, 3, 3, 5, 4] /. 
              symmetricA /. parameters];

Solution:

sol = Solve[Reduce@eqs, allvars][[1]]~Join~parameters~
        Join~{symmetricA}

(AA = doubleA[][{6, 3, 3, diag}, 3] //. sol) // 
      MatrixForm

(SS = doubleS[18, 5] /. sol) // MatrixForm

doubleB[][6, 3, 3, diag] . SS // MatrixForm

Second-order derivative operator:

(DD = der2[18][6, 3, 3, diag] = BuildD2[][{6, 3, 3, diag}, 
              3, 3, 5, 4, sol]) // MatrixForm

Following MN we eliminate all these terms:

ysol = Solve[{DD[[1, 7]] == 0, DD[[1, 8]] == 0, 
            DD[[1, 9]] == 0, DD[[2, 7]] == 0, 
            DD[[2, 8]] == 0, DD[[2, 9]] == 0, 
            DD[[3, 7]] == 0, DD[[3, 8]] == 0, 
            DD[[3, 9]] == 0, DD[[4, 8]] == 0}, {y1, y2, 
            y3, y4, y5, y6, y7, y8, y9, y10}][[1]]

(AAy = AA /. ysol) // MatrixForm

Plot of the eigenvlaues of AA as functions of y1:

eigs = Eigenvalues[AA /. ysol];

Plot[Evaluate[eigs], {y1, - .05,  .05}, 
      PlotRange -> {-1, 1}]

Any negative or small positive value for y1 will do. We take the strange choice in MN (with y1=0 we could decrease the bandwidth of D2).

y1rule = Solve[(DD /. ysol)[[1, 1]] == 114170/40947, 
          y1][[1, 1]]

der2[18][6, 3, 3, diag] = Simplify[DD /. ysol];

Simplify[der2[18][6, 3, 3, diag] /. y1rule] // 
      MatrixForm

checkder2[][6, 3, 3, diag] // MatrixForm

checkS[SS, 18]

Diagonal H, r=8, s=4, tau=4, s2=4, tau3=4, rS=7, tauS=5. Appendix C.4 of MN

First-order operators. The norm has no free parameter:

norm[8, 4, 4, diag] // MatrixForm

The first-derivative has three parameters and this is the (strange) choice of MN:

x123sol = With[{D1 = der1[8, 4, 4, diag]}, 
          Solve[{D1[[1, 2]] == 5544277/5992556, 
              D1[[1, 3]] == 198794991/29962780, 
              D1[[3, 2]] == 9444709/415220}, {x1, x2, 
              x3}]][[1]]

der1[8, 4, 4, diag] /. x123sol // Simplify // MatrixForm

Note that in principle we could get what we want with rS=6. However, for some unknown reason, MN choose rS=7. Free parameters:

parameters = {a[1, 8] -> y1, a[1, 9] -> y2, 
          a[1, 10] -> y3, a[1, 11] -> y4, a[1, 12] -> y5, 
          a[2, 9] -> y6, a[2, 10] -> y7, a[2, 11] -> y8, 
          a[2, 12] -> y9, a[3, 9] -> y10, a[3, 10] -> y11, 
          a[3, 11] -> y12, a[3, 12] -> y13, a[4, 10] -> y14, 
          a[4, 11] -> y15, a[4, 12] -> y16, a[5, 11] -> y17, 
          a[5, 12] -> y18, a[6, 12] -> y19, 
          a[7, 7] -> y20 + 36  y21, a[8, 8] -> y21};

Variables we must solve for:

allvars = 
      Complement[
        Cases[Union[
              Flatten@doubleA[][{8, 4, 4, diag, MN}, 4] /. 
                symmetricA], _a]~Join~
          DeleteCases[Union[Flatten@doubleS[24, 7] /. -x_ -> x], 0], 
        First /@ parameters]

Equations they obey:

eqs = 
        Simplify[D2conditions[][{8, 4, 4, diag, MN}, 4, 4, 7, 5] /. 
              symmetricA /. parameters];

Solution:

sol = Solve[Reduce@eqs, allvars][[1]]~Join~parameters~
        Join~{symmetricA}

(AA = doubleA[][{8, 4, 4, diag, MN}, 4] //. sol) // 
      MatrixForm

This is the (very strange) choice of MN:

S7sol = S[7] -> 6/175;

(SS = doubleS[24, 7] /. sol /. S7sol) // MatrixForm

doubleB[][8, 4, 4, diag] . SS // MatrixForm

Second-order derivative operator:

(DD = der2[24][8, 4, 4, diag, 
              MN] = BuildD2[][{8, 4, 4, diag, MN}, 4, 4, 7, 5, 
              sol~Join~{S7sol}]) // MatrixForm

The corner block has size 8x 12:

DD /. x_Times -> X // MatrixForm

Following MN we eliminate all these coefficients (which do not form a corner):

ysol = Solve[{DD[[1, 8]] == 0, DD[[1, 9]] == 0, 
            DD[[1, 10]] == 0, DD[[1, 11]] == 0, 
            DD[[1, 12]] == 0, DD[[2, 9]] == 0, 
            DD[[2, 10]] == 0, DD[[2, 11]] == 0, 
            DD[[2, 12]] == 0, DD[[3, 9]] == 0, 
            DD[[3, 10]] == 0, DD[[3, 11]] == 0, 
            DD[[3, 12]] == 0, DD[[4, 9]] == 0, 
            DD[[4, 10]] == 0, DD[[4, 11]] == 0, 
            DD[[4, 12]] == 0, DD[[5, 10]] == 0, 
            DD[[5, 11]] == 0, DD[[5, 12]] == 0, 
            DD[[6, 11]] == 0, DD[[6, 12]] == 0, 
            DD[[7, 12]] == 0, DD[[8, 1]] == 0}, {y1, y2, 
            y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14, y15, y16, 
            y17, y18, y19, y20, y21}][[1]]

(AAy = Simplify[AA /. ysol]) // MatrixForm

(DDy = der2[24][8, 4, 4, diag, MN, MN] = 
            Simplify[DD /. ysol]) // MatrixForm

DDy /. x_Times -> X // MatrixForm

We could eliminate one more element from the farthest diagonal taking choosing y21, but not both elements remaining:

Solve[{DDy[[1, 7]] == 0, DDy[[2, 8]] == 0}, {y20, 
        y21}]

Solve[{DDy[[7, 1]] == 0, DDy[[8, 2]] == 0}, {y20, 
        y21}]

The choice in MN:

MNsol = Solve[{DDy[[1, 1]] == 4870382994799/1358976868290, 
            DDy[[2, 2]] == -154646272029/111605472530}, {y20, 
            y21}][[1]]

Note that there is a 0 eigenvalue, as usual:

Eigenvalues[AAy /. MNsol // N] // Chop

I have tried larger A matrices and the fact that all eigenvalues are nonnegative seems robust.

Simplify[DDy /. MNsol] // MatrixForm

checkder2[][8, 4, 4, diag, MN, MN] // MatrixForm

The choice by MN increases the accuracy at one point:

% /. MNsol // MatrixForm

checkS[SS, 24]

Diagonal H, r=8, s=4, tau=4, s2=4, tau3=4, rS=6, tauS=5

First-order operators. The norm has no free parameter:

norm[8, 4, 4, diag] // MatrixForm

The first-derivative has three parameters and this is the (strange) choice of MN:

x123sol = With[{D1 = der1[8, 4, 4, diag]}, 
          Solve[{D1[[1, 2]] == 5544277/5992556, 
              D1[[1, 3]] == 198794991/29962780, 
              D1[[3, 2]] == 9444709/415220}, {x1, x2, 
              x3}]][[1]]

der1[8, 4, 4, diag] /. x123sol // Simplify // MatrixForm

Note that in principle we could get what we want with rS=6. However, for some unknown reason, MN choose rS=7. Free parameters:

parameters = {a[1, 8] -> y1, a[1, 9] -> y2, 
          a[1, 10] -> y3, a[1, 11] -> y4, a[1, 12] -> y5, 
          a[2, 9] -> y6, a[2, 10] -> y7, a[2, 11] -> y8, 
          a[2, 12] -> y9, a[3, 9] -> y10, a[3, 10] -> y11, 
          a[3, 11] -> y12, a[3, 12] -> y13, a[4, 10] -> y14, 
          a[4, 11] -> y15, a[4, 12] -> y16, a[5, 11] -> y17, 
          a[5, 12] -> y18, a[6, 12] -> y19, 
          a[7, 7] -> y20 + 36  y21, a[8, 8] -> y21};

Variables we must solve for:

allvars = 
      Complement[
        Cases[Union[
              Flatten@doubleA[][{8, 4, 4, diag, MN}, 4] /. 
                symmetricA], _a]~Join~
          DeleteCases[Union[Flatten@doubleS[24, 6] /. -x_ -> x], 0], 
        First /@ parameters]

Equations they obey:

eqs = 
        Simplify[D2conditions[][{8, 4, 4, diag, MN}, 4, 4, 6, 5] /. 
              symmetricA /. parameters];

Solution:

sol = Solve[Reduce@eqs, allvars][[1]]~Join~parameters~
        Join~{symmetricA}

(AAbis = doubleA[][{8, 4, 4, diag, MN}, 4] //. sol) // 
      MatrixForm

(SSbis = doubleS[24, 6] /. sol) // MatrixForm

doubleB[][8, 4, 4, diag] . SSbis // MatrixForm

Second-order derivative operator:

(DDbis = der2[24][8, 4, 4, diag, 
              MN] = BuildD2[][{8, 4, 4, diag, MN}, 4, 4, 6, 5, sol]) // 
      MatrixForm

The corner block has size 8x 12:

DDbis /. x_Times -> X // MatrixForm

Following MN we eliminate all these coefficients (which do not form a corner):

ysol = Solve[{DD[[1, 8]] == 0, DD[[1, 9]] == 0, 
            DD[[1, 10]] == 0, DD[[1, 11]] == 0, 
            DD[[1, 12]] == 0, DD[[2, 9]] == 0, 
            DD[[2, 10]] == 0, DD[[2, 11]] == 0, 
            DD[[2, 12]] == 0, DD[[3, 9]] == 0, 
            DD[[3, 10]] == 0, DD[[3, 11]] == 0, 
            DD[[3, 12]] == 0, DD[[4, 9]] == 0, 
            DD[[4, 10]] == 0, DD[[4, 11]] == 0, 
            DD[[4, 12]] == 0, DD[[5, 10]] == 0, 
            DD[[5, 11]] == 0, DD[[5, 12]] == 0, 
            DD[[6, 11]] == 0, DD[[6, 12]] == 0, 
            DD[[7, 12]] == 0, DD[[8, 1]] == 0}, {y1, y2, 
            y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14, y15, y16, 
            y17, y18, y19, y20, y21}][[1]]

(AAybis = Simplify[AAbis /. ysol]) // MatrixForm

(DDybis = der2[24][8, 4, 4, diag, MN, MN] = 
            Simplify[DDbis /. ysol]) // MatrixForm

DDybis /. x_Times -> X // MatrixForm

We could eliminate the whole farthest diagonal, so that the boundary has just one more point than the interior:

ysol2 = Solve[{DDybis[[1, 7]] == 0, 
            DDybis[[2, 8]] == 0, DDybis[[7, 1]] == 0, 
            DDybis[[8, 2]] == 0}, {y20, y21}][[1]]

DDybis /. ysol2 /. x_Times -> X // MatrixForm

Note that there is a 0 eigenvalue, as usual. The minimum bandwidth solution is not valid:

Eigenvalues[AAybis /. ysol2] // N

but there are valid solutions:

Eigenvalues[AAybis /. {y20 -> 0, y21 -> 3}] // N

Simplify[DDy /. {y20 -> 0, y21 -> 3}] // 
      MatrixForm

checkder2[][8, 4, 4, diag, MN, MN] // MatrixForm

The choice by MN increases the accuracy at one point:

% /. {y20 -> 0, y21 -> 3} // MatrixForm

checkS[SS, 24]

Full H, r=4, s=2, tau=3, s2=2, tau2=2, rS=4, tauS=3. Appendix D.1 of MN

First-order operators:

norm[4, 2, 3, MN] // MatrixForm

der1[4, 2, 3, MN] // MatrixForm

This is the operator Q quoted in the paper:

norm[4, 2, 3, MN] . der1[4, 2, 3, MN][[{1, 2, 3, 4}]] // 
        Simplify // MatrixForm

Free parameters:

parameters = {a[1, 1] -> y1, a[1, 2] -> y2, 
          a[1, 3] -> y3, a[1, 4] -> y4, a[1, 5] -> y5, 
          a[1, 6] -> y6, a[2, 1] -> y7, a[2, 2] -> y8, 
          a[2, 3] -> y9, a[2, 4] -> y10, a[2, 5] -> y11, 
          a[2, 6] -> y12, a[3, 5] -> y13, a[3, 6] -> y14, 
          a[4, 5] -> y15, a[4, 6] -> y16};

Variables we must solve for:

allvars = 
      Complement[
        Cases[Union[Flatten@doubleA[][{4, 2, 3}, 2]], _a]~Join~
          DeleteCases[Union[Flatten@doubleS[12, 4] /. -x_ -> x], 0], 
        First /@ parameters]

Equations they obey:

eqs = 
        Simplify[D2conditions[][{4, 2, 3, MN}, 2, 2, 4, 3] /. 
            parameters];

We solve partially for the S variables:

SetAttributes[{a, S}, NHoldAll]

Ssol = Solve[
          eqs[[{1, 2, 3, 4}]], {S[1], S[2], S[3], 
            S[4]}][[1]]

allvars = Complement[allvars, First /@ Ssol]

Still, the system cannot be solved analitically:

sol = 
      ToRules[Reduce[eqs /. Ssol, allvars]]~Join~parameters~Join~
        Ssol

We do it numerically:

sol = 
      Rationalize@ToRules[Reduce[N[eqs /. Ssol, 20], allvars] // Chop]~Join~
        parameters~Join~Ssol

(AA = doubleA[][{4, 2, 3, MN}, 2] //. sol) // 
      MatrixForm

(SS = doubleS[12, 4] /. Ssol) // MatrixForm

(DD = 
          Rationalize[Chop@BuildD2[][{4, 2, 3, MN}, 2, 3, 4, 3, sol]]) // 
      MatrixForm

This is partially the MN case:

DD[[1, 6]]

ysol1 = 
      N[Simplify@Solve[{DD[[1, 5]] == 0, 
                    DD[[1, 6]] == 0} /. Ssol, {y5, y6}][[1]] // 
          ExpandAll, 20]

yeqs = 
      Rationalize[{DD[[1, 1]] == 2, 
                DD[[1, 2]] == -5, DD[[1, 3]] == 4, 
                DD[[1, 4]] == -1} /. Ssol /. ysol1 // 
          Simplify]

ysol2 = 
      Chop@Solve[yeqs[[{1, 2}]], {y11, y12}][[1]]

Chop[yeqs[[{3, 4}]] /. ysol2 // ExpandAll]

ysol3 = Rationalize[
          Solve[%, {y9, y10}] // ExpandAll][[1]]

But there are still other 10 free parameters, that we need to adjust to get A+tr(A) to be a positive matrix.

moreys = 
      Complement[Last /@ parameters, 
        First /@ Join[ysol1, ysol2, ysol3]]

(DDy = DD //. Join[ysol1, ysol2, ysol3] // Simplify // Chop) // 
      MatrixForm

ysol4 = 
      Chop@Solve[{
DDy[[2, 1]] == 
                0.9113401326379418`16, 
DDy[[2, 
                    2]] == -1.639646840154013`16, 
DDy[[2, 3]] == 
                0.4451860342366326`16, 
DDy[[2, 
                    4]] == 
                0.3889216118347602`16, 
DDy[[2, 
                    5]] == -0.1115146289530765`16, 
DDy[[2, 6]] == 
                0.005713690397754591`18, 
DDy[[3, 
                    1]] == -0.02020917311782903`16, 
DDy[[3, 2]] == 
                1.076710314575741`16, 
DDy[[4, 
                    1]] == -0.01972518376035006`16, 
DDy[[4, 2]] == -0.003132697803201818`16}, 
            moreys][[1]]

totalsol = 
      Rationalize[
        Chop@Join[sol /. ysol1 /. ysol2 /. ysol3 /. 
                ysol4, ysol1 /. ysol2 /. ysol3 /. 
                ysol4, ysol2 /. ysol3 /. ysol4, ysol3 /. ysol4, ysol4] // 
          Expand]

(DDyy = der2[12][4, 2, 3, MN] = 
            Rationalize[DDy /. totalsol // Simplify // Chop]) // 
      MatrixForm

(AAyy = Rationalize[AA /. totalsol // Expand]) // 
      MatrixForm

Unfortunately we get negative eigenvalues:

Eigenvalues[AAyy + Transpose[AAyy]] // Chop

Checks:

checkder2[][4, 2, 3, MN] // MatrixForm

checkS[SS, 12]

Full restricted H, r=5, s=2, tau=3, s2=2, tau2=3, rS=5, tauS=3. Peter's case

First-order operators:

norm[5, 2, 3, restr, Peter] // MatrixForm

der1[5, 2, 3, restr, Peter] // MatrixForm

A simple check. This is the operator Q quoted in the paper:

MatrixForm[
      Q = Simplify@
          Rationalize[doublenorm[][5, 2, 3, restr, 
                Peter] . doubleder1[][5, 2, 3, restr, Peter]]]

MatrixForm[Q + Transpose[Q]]

As it is customary, we choose for the second-order operator an interior stencil with the same width as the first-order operator:

MatrixForm[doubleA[][{5, 2, 3, rest, Peter}, 2]]

That leaves us 35 free a-coefficients, plus 5 free S-coefficients. It is not yet clear to me whether we can use just 4 points in S.

We can try to minimize the bandwidth. The optimal situation would drop 13 a-coefficients:

MinBand = 
      Join[Map[# -> 0 &, Apply[a, corner[7, 4], 1]], {a[4, 1] -> 0, 
          a[5, 1] -> 0, a[5, 2] -> 0}]

MatrixForm[doubleA[][{5, 2, 3, rest, Peter}, 2] /. 
        MinBand]

That will turn out to be too restrictive. We leave one point:

MinBand2 = 
      Cases[MinBand, 
        HoldPattern[a[i_, j_] -> 0] /; Abs[i - j] > 3]

MatrixForm[doubleA[][{5, 2, 3, rest, Peter}, 2] /. 
        MinBand2]

There are 29 equations for those 35+4 free coefficients.

eqs = 
        Simplify[D2conditions[][{5, 2, 3, restr, Peter}, 2, 3, 5, 
            3]];

Length[eqs]

We solve partially for four S-variables, leaving S[5] free:

SetAttributes[{a, S}, NHoldAll]

Ssol = Solve[eqs[[{1, 2, 3, 4}]], 
          Array[S, 5]][[1]]

neweqs = 
        DeleteCases[Rationalize@Chop[eqs /. Ssol], True];

and we are left with 25 linear equations for the 35 a-coefficients.

indep = Rationalize@Chop@Reduce[neweqs]

Length[indep]

There are no solutions with a symmetric matrix A:

Reduce[indep /. a[i_, j_] :> a[j, i] /; i > j]

As we said, MinBand is too restrictive. The problem is that equation 16 is inhomogeneous and only involve off-main-diagonals coefficients:

Rationalize@Reduce[List @@ indep /. MinBand]

indep[[16]]

% /. MinBand

%% /. MinBand2

Let us solve for the particular case given by MinBand2:

parameters = {a[3, 6] -> y1, a[4, 6] -> y2, 
          a[4, 7] -> y3, a[5, 7] -> y4};

Rationalize@Reduce[indep /. MinBand2] /. parameters

sol = Join[% // ToRules, parameters, Ssol, MinBand2]

We construct the matrix A with minimum size 10. I have not been able to find a positive case.

MatrixForm[
      AA10 = doubleA[10][{5, 2, 3, rest, Peter}, 2] /. sol /. 
          MinBand2]

This matrix must have nonnegative eigenvalues:

MatrixForm[matrix10 = AA10 + Transpose[AA10]]

It has a nice structure:

matrix10 /. {_Times -> X, _Plus -> X, _Real -> X, 
          y1 -> X} // MatrixForm

minor[range_] := 
      Rationalize@Simplify@
          Det[matrix10[[range, range]]]
minors[size_] := 
      minor /@ Partition[Range[10], size, 1]
allminors = Join @@ Array[minors, 10];

There are solutions considering all but the last minor (54 in total):

Reduce[
      Map[# > 0 &, 
        Drop[allminors, -1] /. {y1 -> 0, y2 -> 0, 
            y3 -> 0}]]

But not with the last minor:

% && (Last[allminors] > 0 /. {y1 -> 0, y2 -> 0, 
              y3 -> 0}) // Reduce

Eigenvalues[
      matrix10 /. {y1 -> 0, y2 -> 0, y3 -> 0, 
          y4 -> 1/4}]

Plot[
      Det[matrix10 /. {y1 -> 0, y2 -> 0, y3 -> 0}], {y4, 0, 
        0.3}]

Let us try the general case:

MatrixForm[
      AA10 = doubleA[10][{5, 2, 3, rest, Peter}, 
            2] /. (indep // ToRules)]

matrix10 = Simplify[AA10 + Transpose[AA10]];

parameters = Union[Cases[AA10, _a | _S, {0, Infinity}]]

minor[range_] := 
      Rationalize@Simplify@
          Det[matrix10[[range, range]]]
minors[size_] := 
      minor /@ Partition[Range[10], size, 1]
allminors = Join @@ Array[minors, 2];

Mathematica cannot even reduce the conditions for the 2-minors...

Reduce[Map[# > 0 &, allminors]] // Simplify // 
      Rationalize

Old stuff:

We construct the operators:

(SS = doubleS[12, 5] //. sol) // MatrixForm

(DD = 
          Rationalize[
            Chop@BuildD2[{5, 2, 3, restr, Peter}, 2, 3, 5, 3, sol] /. 
              example]) // MatrixForm

DD /. {x_Plus -> X, _Times -> X} // MatrixForm

der2[5, 2, 3, restr, Peter] = DD;

Checks:

checkder2[5, 2, 3, restr, Peter] // MatrixForm

Rationalize@checkS[SS, 12]