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

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.

• Maybe this link will be helpful. – Vlad Oct 25 '15 at 9:57
• Are the operators $div$ and ${\nabla}$ taking derivative with respect to $\bf{x}$ or $\bf{y}$? – Hosein Rahnama Oct 28 '15 at 23:21
• with respect to x, i.e. consider that you start from Cartesian coordinates and move to something else – Michael Medvinsky Oct 28 '15 at 23:24
• What is the problem with the final version of @Winther solution? Isn't it satisfactory? – Hosein Rahnama Oct 28 '15 at 23:42
• Can't sleep, because I have the uneasy feeling that what you really want is something completely different. So I could formulate quite another answer tomorrow. Finite Element Method is the keyword. If that is not appropriate either, please let us know. – Han de Bruijn Oct 28 '15 at 23:51

${\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}$$

• I didn't say I'm totally not familiar with tensor calculus, I just said that this didn't help me enough yet. Han's answer doesn't give anything new for me. In all your development A doesn't behave like a matrix, for example $$A\frac{\partial^2 u}{\partial x_i\partial x_i}$$ is a matrix times scalar, but we need it to be scalar. Note that for ($A({\bf x})=\mathrm{diag}(\alpha_1({\bf x}), ...,\alpha_n({\bf x}))$) the translation of variables can be a full matrix.... – Michael Medvinsky Oct 28 '15 at 10:25
• @MichaelMedvinsky Yes I did forget that $A$ was a matrix in your case so I will have to update the answer when I get time. – Winther Oct 28 '15 at 11:16
• see also updated question, with more details and example – Michael Medvinsky Oct 28 '15 at 11:20
• That should fill the bill, I think: (+1) for the competition :-) – Han de Bruijn Oct 28 '15 at 13:39
• I think your expression misses the following term (based on the latest sum) $$a_i G_{\ell i}\frac{\partial G_{k i}}{\partial y_\ell}\frac{\partial u}{\partial y_\ell}$$ What about to write similar expression using scale factors, i.e. $h_{y_i}=\sqrt{\sum_j \left(\frac{\partial x_i}{\partial y_j}\right)^2}$ – Michael Medvinsky Oct 28 '15 at 23:19

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.

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

• yes, polar and spherical coordinates are somewhat not a big hassle, however it doesn't true about elliptic or normal-tangential (local) coordinate system. – Michael Medvinsky Oct 24 '15 at 13:32
• your example is just a Laplacian which is much more simple... – Michael Medvinsky Oct 24 '15 at 13:37
• @MichaelMedvinsky: Unfortunately, I have no experience with the latter two. But the question is: did you actually try tensor calculus? Especially the Einstein convention for repeated indices translates very nicely into loop structures for those indices. – Han de Bruijn Oct 24 '15 at 13:38
• @MichaelMedvinsky: I know that the example given is simple, but I'm promoting an approach here, not the result as such. – Han de Bruijn Oct 24 '15 at 13:41
• the answer is yes and no :) The no is because it is not explicit enough. I can show own work a little, but I'm afraid this will more disturb than help. I've got some book with similar operators, it took me long time to read the Enistein convention and rewrite it to a form that I need and after all it didn't help, may be need another attempt with that book, but this is still in a long term todo list. I would appreciate if you can give in example with the operator I need, may be you can do it better then I've tried. – Michael Medvinsky Oct 24 '15 at 13:48

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.

• 1) you need to use \cdot instead of \dot, 2) what is $\bf A:B$? – Michael Medvinsky Oct 29 '15 at 15:01
• OK, I will add the definition to the answer. :) – Hosein Rahnama Oct 29 '15 at 15:06
• @MichaelMedvinsky: See the answer. :) – Hosein Rahnama Oct 29 '15 at 15:19
• This is nice. I'm not sure, though, that applying the differentiation before the changing of coordinates is the best way. In addition, I don't need derivatives in all directions after changing of variable. May be I wrong, but right now it looks to me redundant work. But still, nice thing. – Michael Medvinsky Oct 29 '15 at 15:55
• @MichaelMedvinsky: Take a look at the edited answer for the final step. :) – Hosein Rahnama Oct 29 '15 at 17:47