The aim is to construct symmetric Laplacian matrix $L$ for 2D square domain.

I have the following discretization of second derivatives on non-uniform grid (will skip the steps of derivation, any textbook works). velocity vectors $u$ and $v$ are defined on cell faces, $x_i$ corresponds to $u_{i+mk}$, where m is amount of unknown $u$'s and k is an integer.

$$ \left(\frac{\partial^2 u}{\partial x^2}\right)_i \approx \frac{u_{i+1}\left(x_i-x_{i-1}\right)-u_i\left(x_{i+1}-x_{i-1}\right)+u_{i-1}\left(x_{i+1}-x_i\right)}{\left(\frac{x_{i+1}-x_{i-1}}{2}\right)\left(x_i-x_{i-1}\right)\left(x_{i+1}-x_i\right)}+O(\Delta x). $$ This scheme becomes second order when $\Delta x_i\approx\Delta x_{i-1}$. The above discretization leads to matrix $\hat{L}$.

Now, consider the symmetrizer matrix $$ \hat{M} =\left[\begin{array}{cc} \frac{1}{2}\left(x_{i+1}+x_{i-1}\right) & 0 \\ 0 & \frac{1}{2}\left(y_{j+1}+ y_{j-1}\right) \end{array}\right], $$ which is a diagonal.

The product $\hat{M}\hat{L}=L$ makes L symmetric for Dirichlet BC and Neumann (1st order) BC.

Now, if I want to use the second order scheme for Neumann BC, discretization becomes 3 point. Therefore, it starts affecting the off diagonal elements in $\hat{L}$ (first oder boundary condition scheme changes only the diagonal element, which does not kill the symmetricity). Hence, $L$ is not symmetric anymore.

Is there a way to modify the diagonal symmetrizer $\hat{M}$ without inducing a lot of computational time? $\hat{M}$ is multiplied by the rhs vector at every time step. I tried the simple $\hat{M}=(\frac{\hat{L}+\hat{L}'}{2})\hat{L}^{-1}$ but it is too dense (see pic below).

M becomes this

Edit: I figured it out. Will post the solution when finish typesetting it in LaTeX.


1 Answer 1


I finally figured it out.

We need to move to the new variables.

(Idk why the typeset result formulas look different than the preview.)

Let us have $m$ intervals in $x$ direction and $n$ intervals in $y$ direction. The number of unknowns is equal to the number of face velocities in the interior of the domain excluding the boundaries. Rewrite $\vec{u}$ as a vector of unknowns

$$\vec{u} = [\underbrace{u_1;u_2;\dotsc;u_{m-1};u_1;u_2;\dotsc;u_{m-1};\dotsc;u_1;u_2;\dotsc;u_{m-1}}_{n\text{ times for each }1...m-1\text{ block}}, \underbrace{v_1;v_2;\dotsc;v_{m};v_1;v_2;\dotsc;v_{m};\dotsc;v_1;v_2;\dotsc;v_{m}}_{n-1 \text{ of each }1,\dotsc,m\text{ block}}].$$

To maintain consistency - indexation goes left to right from bottom to top.

Define matrices

$\Delta y_j=diag([\underbrace{\Delta y_1; \Delta y_1; \dotsc;\Delta y_1}_{m-1\text{ times}}; \underbrace{\Delta y_2; \Delta y_2; \dotsc;\Delta y_2}_{m-1\text{ times}};\dotsc;\underbrace{\Delta y_n; \Delta y_n; \dotsc;\Delta y_n}_{m-1\text{ times}}]),$

$ \Delta x_i=diag([\underbrace{\Delta x_1; \Delta x_2;\dotsc\Delta x_m;\Delta x_1; \Delta x_2;\dotsc\Delta x_m; \dotsc;\Delta x_1; \Delta x_2;\dotsc;\Delta x_m}_{n\text{ times}}])$

\begin{equation} R =\left[\begin{array}{cc} \Delta y_j & 0 \\ 0 & \Delta x_i \end{array}\right], \end{equation}

Let the new vector variable of unknowns (something like mass flux) be defined as

$$\vec{q} = R\vec{u},$$

where $R$ is a diagonal matrix of size $(m-1)n\times m(n-1)$ defined above.

Next, we need to define the difference matrices to cancel out the $\frac{x_{i+1}−x_{i-1}}{2}$ central term in diffusion discretization.


${\Delta x_i+\Delta x_{i-1}} = diag[\underbrace{\Delta x_2+\Delta x_1; \Delta x_3+\Delta x_2; \dotsc \Delta x_m-\Delta x_{m-1};...;\Delta x_2+\Delta x_1; \Delta x_3+\Delta x_2; \dotsc \Delta x_m-\Delta x_{m-1}}_{n\text{ times for each block of } \Delta x_2+\Delta x_1\text{ to }\Delta x_m-\Delta x_{m-1}}]$

${\Delta y_j+\Delta y_{j-1}} = diag[\underbrace{\Delta y_2+\Delta y_1; \dotsc \Delta y_2+\Delta y_1;}_{m\text{ times}} \underbrace{\Delta y_3+\Delta y_2; \dotsc \Delta y_3+\Delta y_2;}_{m\text{ times}} \underbrace{\Delta y_n+\Delta y_{n-1}\dotsc \Delta y_n+\Delta y_{n-1}}_{m\text{ times}}],$

\begin{equation} \hat{M} =\left[\begin{array}{cc} \frac{1}{2}\left(\Delta x_i+\Delta x_{i-1}\right) & 0 \\ 0 & \frac{1}{2}\left(\Delta y_j+\Delta y_{j-1}\right) \end{array}\right]. \end{equation}

Now, since $\vec{q}=R\vec{u}$ then $\vec{u}=R^{-1}\vec{q}$, Going back to Laplacian matrix

$$L\vec{u}= LR^{-1}\vec{q}.$$

Premultiply by $\hat M$, the expression will become

$$\hat MLR^{-1}\vec{q},$$

where $\hat M L R^{-1}$ is symmetric by construction.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.