Cómo simplificar $\frac{\partial^m}{\partial y_i^m}\mathrm{div }(A\nabla u({\bf x}({\bf y})))$

Question

Dejemos que $\bf{x}\in\mathbb{R}^n$ (interesante en $n\in\{2,3\}$ ) y que $A=A_{n\times n}=\mathrm{diag}(a_1({\bf x}),\dots,a_n({\bf x}))$ Es decir $$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} $$ Supongamos que $a_i\ne a_j, \forall 1\le i,j\le n$ .

Dejemos que $\bf y$ sean algunas coordenadas curvilíneas, por ejemplo, elípticas, esféricas, normal-tangenciales. Me interesa diferenciar $$\mathrm{L}u({\bf x}) =\{\mathrm{div }(A\nabla u)\}({\bf x})$$ con respecto a $y_i$ es decir, expresar $$\frac{\partial^m}{\partial y_i^m} \mathrm{L}u({\bf x}),\quad \text{interesting in no more then }0\le m\le4$$ de forma explícita pero compacta (por $m=0$ Me refiero a un simple cambio de coordenadas).

La expresión, que obtuve hasta ahora (usando factores de escala \metrics ), se vuelven demasiado largas y feas. Digamos que se vuelve bastante inmanejable cuando necesito programarlo, lo que lo hace vulnerable para los bugs y puede afectar la estabilidad numérica también.

Ejemplo en 2D

Dejemos que $a_1=a, a_2=b$ uno traduce $$ \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} $$ a unas coordenadas (s,t) para obtener $$ \alpha u_s + \beta u_t +\gamma u_{st} +\delta u_{tt} +\sigma u_{ss} $$ y diferenciarlo, por un instante con respecto a $s$ para conseguir $$\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}$$

Se agradecerá cualquier ayuda.

Answer 1

${\bf NB:}$ Esta respuesta es para la pregunta original antes de ser editada: cómo cambiar las coordenadas de la expresión $\text{div}(A\nabla u)$ y ponerlo de forma explícita y compacta.

---

Aquí hay una derivación muy básica de esa forma explícita y relativamente compacta que pides y que permite una implementación numérica bastante sencilla. Existen formas mejores (más sistemáticas y estructuradas) de hacerlo (véase la respuesta de Han), pero el objetivo de esta respuesta es hacerlo de la forma más básica posible y evitar introducir demasiadas novedades. No obstante, utilizaré La convención de la suma de Einstein : cada uno de los índices repetidos debe ser sumado sobre $1,2,\ldots,n$ en las derivaciones siguientes.

---

Bajo un cambio de variables de ${\bf x} = \{x_1,\ldots,x_n\} \to {\bf y} = \{y_1,\ldots,y_n\}$ las derivadas parciales se transforman en $$\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}$$

donde he definido la matriz $G_{ji} \equiv \frac{\partial y_j}{\partial x_i}$ . Tenemos además que las segundas derivadas se transforman como

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

La cantidad que se busca se puede escribir

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

que se expresa en ${\bf y}$ coordenadas se convierte en

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

o con una suma explícita

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

Esta forma parece bastante sencilla de aplicar numéricamente. Para cada sistema de coordenadas se pueden calcular las componentes de $G_{ji}$ y sus derivadas de forma analítica y almacenarlos en una función que tome $i,j$ y ${\bf y}$ como argumentos. La página web $G_{ji}$ suelen ser expresiones bastante sencillas. No olvide que ${\bf x}$ debe expresarse en ${\bf y}$ que entra en los argumentos de $A$ en la fórmula anterior.

Desde su $A$ la matriz es diagonal $A = \text{diag}(a_1,\ldots,a_n)$ podemos simplificar aún más

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

Es posible simplificar aún más con algo de trabajo (pero no sé hasta qué punto será útil). Si dejamos que $H$ sea el Hessian de $u$ entonces podemos, por ejemplo, escribir en términos de operaciones matriciales lo siguiente

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

Los otros términos son un poco más complicados ya que implican derivadas de una matriz (estos términos corresponden a Símbolos de Christoffel ), por lo que no se puede representar en la forma de matriz 2D estándar.

---

Es muy fácil cometer pequeños errores en este tipo de cálculos, así que vamos a hacer una sencilla prueba para $n=2$ y coordenadas polares para ver que funciona al menos para este caso. Tenemos ${\bf x} = (x,y) = (r\cos\theta,r\sin\theta)$ y ${\bf y} = (r,\theta)$ así que

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

Evaluando la fórmula anterior (para $a_i = 1$ así que $A$ es la identidad) utilizando 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]

nos da el resultado correcto

$$\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 2

A primera vista parece algo descabellado, pero una vez que se dominen algunos primeros principios, se verá que Cálculo tensorial puede ser indispensable para que estas cosas sean susceptibles de un tratamiento sistemático. No es fácil empezar, pero mi experiencia personal es que vale la pena.
Aquí hay algunas referencias de Wikipedia, pero puede ser preferible un buen libro de texto.

• Cálculo de Ricci
• Tensores en coordenadas curvilíneas

A continuación se muestra una derivación del gradiente y el operador de Laplace, para coordenadas esféricas, como se encuentra en mis viejos apuntes de la universidad:

Answer 3

Represento mi solución en algunos pasos. Como puede resultar demasiado largo, dejo la demostración de algunas partes para ti. Además, en lo que sigue, $u$ es un escalar, ${\bf{v}}$ , ${\bf{a}}$ , ${\bf{b}}$ son vectores, ${\bf{A}}$ , ${\bf{B}}$ son tensores (o matrices) de segundo orden. Además, hay que tener en cuenta las siguientes definiciones

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

donde la suma está implícita sobre los índices repetidos. Estas definiciones se pueden encontrar en cualquier libro de cálculo tensorial estándar. De hecho, deberías estar familiarizado con $:$ el producto escalar de dos tensores de segundo orden y con $\otimes $ El producto tensorial, y se pueden derivar las fórmulas anteriores por cálculo directo.

Paso 1 . La siguiente fórmula es válida

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

En esta fórmula y en lo que sigue, se puede poner ${\bf{v}} = \nabla u$ para conseguir lo que quieres.

Paso 2 . Considere estas fórmulas

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

He demostrado estas fórmulas por inducción. Son algunas generalizaciones de la conocida fórmula

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

Paso 3 . Utilice $(1)$ y $(2)$ para obviar

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

y algunas simplificaciones claras llevarán a

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

Paso 4 . El último paso es aplicar la regla de la cadena. En todas las fórmulas anteriores se ha supuesto que $\nabla = {\nabla _{\bf{x}}}$ . Se puede demostrar que

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

Sin embargo, no tuve tiempo de ir más allá. Pero debería derivar algunas fórmulas basadas en $(5)$ para obtener algo como

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

y luego combinar $(6)$ con $(4)$ para obtener el resultado final. Pero esto no será fácil. Para obtener más ayuda sobre este último paso, consulte este enlace .