4

I am working to reproduce results from a paper [links directly to page 17, equations 3.4 on some browsers] on a finite-volume method reconstruction method called WENO which I'm using in one dimension. Part of the method involves calculating a so-called smoothness indicator $\beta$ for cell average values $\bar{u}_i$ in a stencil which consists of a small number of cells, in my case three. For uniform meshes the paper gives explicit formulae for $\beta$ in terms of the cell average values. For example, a result they give for a stencil three cells wide is,

$$\beta_1 = \frac{13}{12}\left(\bar{u}_{i - 2} - 2\bar{u}_{i - 1} + \bar{u}_i\right)^2 + \frac{1}{4}\left(\bar{u}_{i-2} - 4\bar{u}_{i - 1} + 3\bar{u}_i\right)^2$$

The other formulae they give for various shifts of stencils of three cells all have this same form, albeit with different coefficients for the cell average terms within the parentheses. [Does this form have a name? If so, I'd like to improve the question title.]

The software I have written, which uses the computer algebra system Sympy to perform the necessary calculus to get to this point, produces the following expression which is equivalent to the above, and so correct,

$$\beta_1 = \frac{10}{3}{\bar{u}_i}^2 - \frac{31}{3}\bar{u}_{i}\bar{u}_{i-1} + \frac{11}{3}\bar{u}_{i}\bar{u}_{i-2} + \frac{25}{3}{\bar{u}_{i-1}}^2 - \frac{19}{3}\bar{u}_{i-1}\bar{u}_{i-2} + \frac{4}{3}{\bar{u}_{i-2}}^2$$

What method or steps can I follow to rearrange my result to the one given in the paper? Ideally I'm looking for a method with generalises to the other equations in the paper which have a similar structure.

Rob Smallshire
  • 93
  • 7
  • Your link to the paper doesn't show the relevant page, but makes to download the full paper. Could you provide a link to the page? – Viera Čerňanová Sep 09 '21 at 21:30
  • Is there a typo in your second equation? [e.g. should the +s in $\bar{u}_{i+1}$ be -s instead?] – John Joy Sep 10 '21 at 00:29
  • @JohnJoy Yes. My software produces two equations which are identical except for these signs and I mistakenly used the wrong one. Thank you. I've updated the question. – Rob Smallshire Sep 10 '21 at 06:44
  • 1
    @user376343 My browser understands the page link, yours apparently does not. Sorry for not being more explicit. The relevant equations are at the top of page 17 in the paper, in equations 3.4. – Rob Smallshire Sep 10 '21 at 06:45
  • A related question is to know why the authors chose this particular form. – Rob Smallshire Sep 12 '21 at 09:21
  • 1
    The name you're looking for might be “sum of squares” , or in this case weighted sum of squares. Its value is in guaranteeing that output values are always non-negative (assuming in the weighted case that the weights are all non-negative, which in yr paper they are). It might have no further importance beyond that? The hint I see is bottom of p16 “The smoothness indicator βl defined in (3.3) is a scaled square sum of the L2- norms …”.

    If so, your next step might to look for a Python Sum Of Squares Solver ?

     – Chris F Carroll Sep 12 '21 at 15:49
  • 1
    One possible way to start is to look up algorithms for how to write polynomials in depressed form. For instance, if you use the Mathematica code here, then, for instance, Simplify /@ depressPolynomial[ 10/3 x^2 - 31/3 x y + 11/3 x z + 25/3 y^2 - 19/3 y z + 4/3 z^2, z] yields $\frac{1}{48} (11 x-19 y+8 z)^2+\frac{13}{16} (x-y)^2$. Using z or y instead yields a different sum of squares. So the form is not unique, but it might get you started. – march Sep 12 '21 at 16:13
  • I'm not an expert but you might like to try simplify() from Sympy https://docs.sympy.org/latest/tutorial/simplification.html . The problem is that computer algebra systems are not magic, nor human intelligent, and essentially try a set of heuristics to simplify. If it works, great, if not you might have to contact the development team of sympy at GitHub via their issues to see if this is a case they can implement. Mathematica is far more mature (and expensive) and would likely to do this simplification. -- ultimately, if the expressions are equivalent, do you need them simplified? – Penelope Sep 12 '21 at 22:07
  • 1
    @Tariq. Sympy can't do what I want out of the box. I think I might need a sum of squares factorisation algorithm, as mentioned here with respect to Sympy: https://github.com/sympy/sympy/issues/11332 – Rob Smallshire Sep 13 '21 at 17:39

1 Answers1

1

In the case that the weights are all non-negative, the coefficient matrix $A$ of the quadratic form is positive semi-definite. For instance, if we consider $3\beta_1$ in your case, it's easy to check $$ A = \begin{pmatrix} 10 & -31/2 & 11/2 \\ -31/2 & 25 & -9/2 \\ 11/2 & -9/2 & 4 \\ \end{pmatrix} $$ has eigenvalues of $39\pm3\sqrt{39}$ and $0$, thus it's positive semi-definite.

There are numerous ways to get the decomposition (or square root) of $A$, i.e., if one has $A=B^TB$, then $B\bar{u}$ gives each square term.

I don't know Sympy, but one famous method is Cholesky decomposition, which Sympy likely has an implementation. Below is one such decomposition (by using Wolfram Alpha): $$3\beta_1=\frac{1}{40}(20\bar{u}_i-31\bar{u}_{i-1}+11\bar{u}_{i-2})^2+\frac{39}{40}(\bar{u}_{i-1}-\bar{u}_{i-2})^2$$

Wiley
  • 1,437