# 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.

Oct 25, 2015 at 9:57
• Are the operators $div$ and ${\nabla}$ taking derivative with respect to $\bf{x}$ or $\bf{y}$? Oct 28, 2015 at 23:21
• with respect to x, i.e. consider that you start from Cartesian coordinates and move to something else Oct 28, 2015 at 23:24
• What is the problem with the final version of @Winther solution? Isn't it satisfactory? Oct 28, 2015 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. Oct 28, 2015 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.... Oct 28, 2015 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. Oct 28, 2015 at 11:16
• see also updated question, with more details and example Oct 28, 2015 at 11:20
• That should fill the bill, I think: (+1) for the competition :-) Oct 28, 2015 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}$ Oct 28, 2015 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. Oct 24, 2015 at 13:32
• your example is just a Laplacian which is much more simple... Oct 24, 2015 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. Oct 24, 2015 at 13:38
• @MichaelMedvinsky: I know that the example given is simple, but I'm promoting an approach here, not the result as such. Oct 24, 2015 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. Oct 24, 2015 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$? Oct 29, 2015 at 15:01
• OK, I will add the definition to the answer. :) Oct 29, 2015 at 15:06
• @MichaelMedvinsky: See the answer. :) Oct 29, 2015 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. Oct 29, 2015 at 15:55
• @MichaelMedvinsky: Take a look at the edited answer for the final step. :) Oct 29, 2015 at 17:47