How to simplify $\frac{\partial^m}{\partial y_i^m}\mathrm{div }(A\nabla u({\bf x}({\bf y})))$

Question

Let $\bf{x}\in\mathbb{R}^n$ (interesting in $n\in\{2,3\}$) and let $A=A_{n\times n}=\mathrm{diag}(a_1({\bf x}),\dots,a_n({\bf x}))$, that is $$A_{2\times2}=\begin{bmatrix} a_1(\bf{x})&0\\0&a_2(\bf{x}) \end{bmatrix} \quad \text{or}\quad A_{3\times3}=\begin{bmatrix} a_1(\bf{x})&0&0\\0&a_2(\bf{x})&0\\ 0&0& a_3(\bf{x}) \end{bmatrix} $$ Assume that $a_i\ne a_j, \forall 1\le i,j\le n$.

Let $\bf y$ be some curvilinear coordinates, e.g. elliptic, spherical, normal-tangential. I'm interesting to differentiate $$\mathrm{L}u({\bf x}) =\{\mathrm{div }(A\nabla u)\}({\bf x})$$ with respect to $y_i$ , that is to express $$\frac{\partial^m}{\partial y_i^m} \mathrm{L}u({\bf x}),\quad \text{interesting in no more then }0\le m\le4$$ explicitly but compact (by $m=0$ I mean just change of coordinates).

The expression, that I got so far (using scale factors\metrics), become too long and ugly. Say it is become pretty much unmanageable when I need to program it, which make it vulnerable for bugs and may affect numerical stability as well.

Example in 2D

Let $a_1=a, a_2=b$ one translates $$ \mathrm{L} u = \frac{\partial}{\partial x} \left( a u_x \right) + \frac{\partial}{\partial y} \left( b u_y \right) = a_x u_x + a u_{xx} + b_y u_y + b u_{yy} $$ to some (s,t) coordinates to get $$ \alpha u_s + \beta u_t +\gamma u_{st} +\delta u_{tt} +\sigma u_{ss} $$ and differentiate it, for an instant with respect to $s$ it to get $$\begin{align} & \alpha_s u_s +\alpha u_{ss} + \beta_s u_t + \beta u_{ts} +\gamma_s u_{st}+\gamma u_{sst} +\delta_s u_{tt}+\delta u_{tts} +\sigma_s u_{ss}+\sigma u_{sss}\\ =& \alpha_s u_s + \beta_s u_t+ (\alpha +\sigma_s) u_{ss} +\delta_s u_{tt}+ (\beta +\gamma_s) u_{st} +\gamma u_{sst}+\delta u_{tts}+\sigma u_{sss} \end{align}$$

Any help will be appreciated.

Answer 1

It seems to be somewhat far-fetched at first sight, but once you master some first principles, it will become apparent that Tensor Calculus may be indispensable to make these things amenable to systematic treatment. It's not easy to start with, but my personal experience is that it pays off.
Here are some Wikipedia references, but a good textbook may be preferable.

• Ricci calculus
• Tensors in curvilinear coordinates

Below is a derivation of the gradient & Laplace operator, for spherical coordinates, as found in my old college notes:

Answer 2

${\bf NB:}$ This answer is for the original question before it was edited: how to change coordinates of the expression $\text{div}(A\nabla u)$ and put it in an explicit and compact form.

---

Here is a very basic derivation of such a explicit and relatively compact form that you ask for which allows for a fairly simple implementation numerically. There do exist better (more systematic and structured) ways to do this (see Han's answer), but the aim of this answer is to just do it as basic as possible and avoid introducing too much new. I will however use Einsteins summation convention: every repeated indices are to be summed over $1,2,\ldots,n$ in the derivations below.

---

Under a change of variables from ${\bf x} = \{x_1,\ldots,x_n\} \to {\bf y} = \{y_1,\ldots,y_n\}$ the partial derivatives transform as $$\frac{\partial }{\partial x_i} = \frac{\partial y_j}{\partial x_i}\frac{\partial}{\partial y_j} \equiv G_{ji}\frac{\partial}{\partial y_j}$$

where I have defined the matrix $G_{ji} \equiv \frac{\partial y_j}{\partial x_i}$. We further have that the second derivatives transform as

$$\frac{\partial^2 }{\partial x_i\partial x_{j}} = G_{\ell i}G_{kj}\frac{\partial^2}{\partial y_\ell\partial y_k} + \frac{\partial G_{\ell i}}{\partial x_{\ell}}\frac{\partial}{\partial y_\ell} = G_{\ell i}G_{kj}\frac{\partial^2}{\partial y_\ell\partial y_k} + G_{kj}\frac{\partial G_{\ell i}}{\partial y_k}\frac{\partial}{\partial y_\ell}$$

The quantity you are after can be written

$$\text{div}(A\nabla u) \equiv \frac{\partial}{\partial x_i}\left(A_{ij}\frac{\partial u}{\partial x_j}\right) = A_{ij}\frac{\partial^2 u}{\partial x_i\partial x_j} + \frac{\partial A_{ij}}{\partial x_i}\frac{\partial u}{\partial x_j}$$

which expressed in ${\bf y}$ coordinates becomes

$$\text{div}(A\nabla u) = A_{ij}\left[G_{\ell i}G_{kj}\frac{\partial^2 u}{\partial y_\ell\partial y_k} + G_{kj}\frac{\partial G_{\ell i}}{\partial y_k}\frac{\partial u}{\partial y_\ell}\right] + G_{\ell i}G_{k j}\frac{\partial A_{ij}}{\partial y_\ell}\frac{\partial u}{\partial y_k}$$

or with explicit summation

$$\text{div}(A\nabla u)({\bf y}) = \sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^n\sum_{\ell=1}^n\left\{ A_{ij}\left[G_{\ell i}G_{kj}\frac{\partial^2 u}{\partial y_\ell\partial y_k} + G_{kj}\frac{\partial G_{\ell i}}{\partial y_k}\frac{\partial u}{\partial y_\ell}\right] + G_{\ell i}G_{k j}\frac{\partial A_{ij}}{\partial y_\ell}\frac{\partial u}{\partial y_k} \right\}$$

This form looks fairly simple to implement numerically. For each coordinate-system you can compute the components of $G_{ji}$ and it's derivatives analytically and store them in a function that takes $i,j$ and ${\bf y}$ as arguments. The $G_{ji}$'s will usually be fairly simple expressions. Don't forget that ${\bf x}$ should be expressed in ${\bf y}$ coordinates which goes into the arguments of $A$ in the formula above.

Since your $A$ matrix is diagonal $A = \text{diag}(a_1,\ldots,a_n)$ we can simplify even further

$$\text{div}(A\nabla u)({\bf y}) = \sum_{i=1}^n\sum_{k=1}^n\sum_{\ell=1}^n\left\{ a_i\left[G_{\ell i}G_{ki}\frac{\partial^2 u}{\partial y_\ell\partial y_k} + G_{ki}\frac{\partial G_{\ell i}}{\partial y_k}\frac{\partial u}{\partial y_\ell}\right] + G_{\ell i}G_{k i}\frac{\partial a_i}{\partial y_\ell}\frac{\partial u}{\partial y_k} \right\}$$

Even further simplifications are possible with some work (but I don't know how useful it will be). If we let $H$ be the Hessian of $u$ then we can for example write in terms of matrix operations as follows

$$ A_{ij}G_{\ell i}G_{kj}\frac{\partial^2 u}{\partial y_\ell\partial y_k} = \text{Tr}[GAG^TH]$$

The other terms are a bit more tricky as they involve derivatives of a matrix (these terms corresponds to Christoffel symbols) making it not possible to represent in the standard 2D matrix form.

---

It's really easy to make small mistakes in such calculations so let's do a simple test for $n=2$ and polar coordinates to see that it works at least for this case. We have ${\bf x} = (x,y) = (r\cos\theta,r\sin\theta)$ and ${\bf y} = (r,\theta)$ so

$$\pmatrix{&G_{11} = \frac{\partial r}{\partial x} = \cos\theta &G_{12} = \frac{\partial r}{\partial y} = \sin\theta\\ &G_{21} = \frac{\partial \theta}{\partial x} = -\frac{\sin\theta}{r} &G_{22} = \frac{\partial \theta}{\partial y} = \frac{\cos(\theta)}{r} } $$

Evaluating the formula above (for $a_i = 1$ so $A$ is the identity) using mathematica

coord = {r, t};
g = {{Cos[t], Sin[t]}, {-Sin[t]/r, Cos[t]/r}};
a = {1, 1, 1};
laplacian = 0;
Do[
  laplacian += 
     a[[i]] (  g[[l, i]] g[[k, i]] D[u[r, t], coord[[l]], coord[[k]]] 
            +  g[[k, i]] D[g[[l, i]], coord[[k]]] D[u[r, t], coord[[l]]] ) 
            + g[[l, i]] g[[k, i]] D[a[[i]], coord[[l]]] D[u[r, t], coord[[k]]];
, {k, 1, 2}, {l, 1, 2}, {i, 1, 2}]
FullSimplify[laplacian]

gives us the correct result

$$\nabla^2 u = \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial \theta^2}$$

Answer 3

I represent my solution in some steps. As it may become too lengthy, I leave the proof of some parts to you. Also, in what follows, $u$ is a scalar, ${\bf{v}}$, ${\bf{a}}$, ${\bf{b}}$ are vectors, ${\bf{A}}$, ${\bf{B}}$ are second order tensors (or Matrices). Furthermore, you should notice the following definitions

$$\begin{array}{l} {\bf{A}}:{\bf{B}} = {A_{ij}}{B_{ij}}\\ \nabla .{\bf{A}} = \frac{{\partial {A_{ij}}}}{{\partial {x_i}}}{{\bf{e}}_j} \end{array}$$

where summation is implied over repeated indices. These definitions can be found in any standard tensor calculus book. In fact, you should be familiar with $:$, the scalar product of two second order tensors and with $\otimes $, the tensor product, and you can just derive the above formulas by direct computation.

Step 1. The following formula is valid

$$\begin{array}{l} \nabla .\left( {{\bf{A}}.{\bf{v}}} \right) = \left( {\nabla .{\bf{A}}} \right).{\bf{v}} + {\bf{A}}:\nabla {\bf{v}}\\\end{array}\tag{1}$$

In this formula and what follows, you can put ${\bf{v}} = \nabla u$ to get what you want.

Step 2. Consider these formulas

$$\begin{array}{l} {\nabla ^{\left( n \right)}}\left( {{\bf{a}}.{\bf{b}}} \right) = \sum\limits_{k = 0}^{n - 1} {\left( {\begin{array}{*{20}{c}} {n - 1}\\ k \end{array}} \right)\left( {{\nabla ^{\left( {n - k} \right)}}{\bf{a}}.{\nabla ^{\left( k \right)}}{\bf{b}} + {\nabla ^{\left( {n - k} \right)}}{\bf{b}}.{\nabla ^{\left( k \right)}}{\bf{a}}} \right)} ,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,n \ge 1\\ {\nabla ^{\left( n \right)}}\left( {{\bf{A}}:{\bf{B}}} \right) = \sum\limits_{k = 0}^{n - 1} {\left( {\begin{array}{*{20}{c}} {n - 1}\\ k \end{array}} \right)\left( {{\nabla ^{\left( {n - k} \right)}}{\bf{A}}:{\nabla ^{\left( k \right)}}{\bf{B}} + {\nabla ^{\left( {n - k} \right)}}{\bf{B}}:{\nabla ^{\left( k \right)}}{\bf{A}}} \right)} ,\,\,\,\,\,\,\,\,n \ge 1\\ {\nabla ^{\left( n \right)}} = \underbrace {\nabla \nabla ...\nabla }_n \end{array}\tag{2}$$

I proved theses formulas by induction. These are some generalization of the well known formula

$${\left( {fg} \right)^{(n)}} = \sum\limits_{k = 0}^n {\left( {\begin{array}{*{20}{c}} n\\ k \end{array}} \right){f^{(n - k)}}{g^{(k)}}}$$

Step 3. Use $(1)$ and $(2)$ to obtian

$$\begin{array}{l} {\nabla ^{\left( n \right)}}\left[ {\nabla .\left( {{\bf{A}}.{\bf{v}}} \right)} \right] = {\nabla ^{\left( n \right)}}\left[ {\left( {\nabla .{\bf{A}}} \right).{\bf{v}}} \right] + {\nabla ^{\left( n \right)}}\left[ {{\bf{A}}:\nabla {\bf{v}}} \right]\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \sum\limits_{k = 0}^{n - 1} {\left( {\begin{array}{*{20}{c}} {n - 1}\\ k \end{array}} \right)\left( {{\nabla ^{\left( {n - k} \right)}}\left( {\nabla .{\bf{A}}} \right).{\nabla ^{\left( k \right)}}{\bf{v}} + {\nabla ^{\left( {n - k} \right)}}{\bf{v}}.{\nabla ^{\left( k \right)}}\left( {\nabla .{\bf{A}}} \right)} \right)} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + \sum\limits_{k = 0}^{n - 1} {\left( {\begin{array}{*{20}{c}} {n - 1}\\ k \end{array}} \right)\left( {{\nabla ^{\left( {n - k} \right)}}{\bf{A}}:{\nabla ^{\left( k \right)}}\nabla {\bf{v}} + {\nabla ^{\left( {n - k} \right)}}\nabla {\bf{v}}:{\nabla ^{\left( k \right)}}{\bf{A}}} \right)} \end{array}\tag{3}$$

and some clear simplifications will lead to

$${\nabla ^{\left( n \right)}}\left[ {\nabla .\left( {{\bf{A}}.{\bf{v}}} \right)} \right] = \sum\limits_{k = 0}^{n - 1} {\left( {\begin{array}{*{20}{c}} {n - 1}\\ k \end{array}} \right)\left( \begin{array}{l} {\nabla ^{\left( {n - k} \right)}}\left( {\nabla .{\bf{A}}} \right).{\nabla ^{\left( k \right)}}{\bf{v}} + {\nabla ^{\left( {n - k} \right)}}{\bf{v}}.{\nabla ^{\left( k \right)}}\left( {\nabla .{\bf{A}}} \right)\\ + {\nabla ^{\left( {n - k} \right)}}{\bf{A}}:{\nabla ^{\left( {k + 1} \right)}}{\bf{v}} + {\nabla ^{\left( {n - k + 1} \right)}}{\bf{v}}:{\nabla ^{\left( k \right)}}{\bf{A}} \end{array} \right)}\tag{4}$$

Step 4. The last step is to apply the chain-rule. In all above formulas it was assumed that $\nabla = {\nabla _{\bf{x}}}$. It can be proved that

$$\begin{array}{*{20}{l}} {{\nabla _{\bf{y}}}{\bf{v}} = {\nabla _{\bf{y}}}{\bf{x}}.{\nabla _{\bf{x}}}{\bf{v}}}\\ {{\nabla _{\bf{y}}}{\bf{A}} = {\nabla _{\bf{y}}}{\bf{x}}.{\nabla _{\bf{x}}}{\bf{A}}} \end{array}\tag{5}$$

However, I didn't have time to go further. But you should derive some formulas based on $(5)$ to obtain something like

$$\nabla _{\bf{y}}^{(n)}\left[ {{\nabla _{\bf{y}}}.\left( {{\bf{A}}.{\bf{v}}} \right)} \right] = \left( {..?..} \right)\nabla _{\bf{x}}^{(n)}\left[ {{\nabla _{\bf{x}}}.\left( {{\bf{A}}.{\bf{v}}} \right)} \right]\tag{6}$$

and then combine $(6)$ with $(4)$ to derive the final result. But this won't be easy! For further help on this last step see this link.