Finite difference approximation of elliptic operator: divergence form vs. non-divergence form

536 Views Asked by At

Consider the principal part of a standard, strongly elliptic linear operator of 2nd order. In divergence form, this is given by $$ \nabla \cdot (a \nabla u), $$ while in non-divergence form, it is given by expanding out this expression: $$ a \nabla^2 u + \nabla a \cdot \nabla u $$ where both $a,u$ are smooth, scalar functions, and $a$ is always positive. We seek to discretize this operator using finite differences on a regular grid. Replacing $\nabla$ with $D$, the finite difference operator, the question is, when do we want to use the divergence form, versus the non-divergence form. In one dimension, for example, applying central differences to the divergence form yields $$ h^{-2} \left( a_{1/2} (u_1 - u_0) - a_{-1/2} (u_0 - u_{-1}) \right). $$ Meanwhile, applying the same to the non-divergence form yields $$ h^{-2} \left(a_0 (u_1 - 2u_0 + u_{-1}) + (a_1 - a_{-1})(u_1 - u_{-1}) \right). $$ It has been impressed upon me by a colleague that the divergence form is preferable because the discretization gives rise to a symmetric matrix, which matches the self-adjoint-ness of the original differential operator. From the standpoint of linear solvers, it is clear that this is advantageous. Indeed, for Galerkin formulations, it is clear why you want to retain the divergence form since then you can provide a good weak formulation, and this has been how I've approached this problem in the past.

However, if we ignore the fact that symmetric linear solvers have better performance, it is not clear to me what numerical advantages the divergence form has. Indeed, both methods are second-order, and if you explicitly expand out the truncation error, the non-divergence form has smaller truncation error term, since the chain rule expansion was done analytically, leading to fewer error terms.

From an approximation theory point of view, what exactly makes the divergence form preferable (if at all)? How do we formalize the intuition that preserving symmetricity in the discretization matrix will lead to a better numerical approximation?