PDEs elípticas

You are currently browsing articles tagged PDEs elípticas.

En el libro de Análisis numérico de Burden y Faires aparecen una serie de ejemplos y ejercicios resueltos de PDEs elípticas en 2D. Vamos a intentar resolverlos utilizando Chombo…

La primera ecuación corresponde con un ejemplo y es

$latex frac{partial^2}{partial x^2}u(x,y) + frac{partial^2}{partial y^2}u(x,y) = 0$

en

$latex Omega = { (x,y): 0<x<0.5, 0<y<0.5}$

y cuyas condiciones en la frontera $latex partial Omega$ son

$latex u(x,0)=0, u(0,y)=0, u(x,0.5)=200x, u(0.5,y)=200y$.

El resultado es:

ejeM1

Para el ejercicio 12.1.1

$latex u_{xx} + u_{yy} = 4$

en

$latex Omega = { (x,y): 0<x<1, 0<y<2}$

y cuyas condiciones en la frontera $latex partial Omega$ son

$latex u(x,0)=x^2, u(0,y)=y^2, u(x,2)=(x-2)^2, u(1,y)=(y-1)^2$.

El resultado es:

ejeR1

A continuación, para el 12.1.3.a

$latex u_{xx} + u_{yy} = 0$

en

$latex Omega = { (x,y): 0<x<1, 0<y<1}$

y cuyas condiciones en la frontera $latex partial Omega$ son

$latex u(x,0)=0, u(0,y)=0, u(x,1)=x, u(1,y)=y$.

El resultado es:

ejeR3a

En el 12.1.3.b encontramos

$latex Delta u = -(cos(x+y)+cos(x-y))$

en

$latex Omega = { (x,y): 0<x<pi, 0<y<frac{pi}{2}}$

y

$latex u(x,0)=cos x , u(0,y)=cos y, u(x,frac{pi}{2})=0, u(pi,y)=-cos y$,

obteniendo:

 ejeR3b ejeR3b2

Tags: , ,

En $latex 1D$ tenemos:

$latex u_{xx} = f(x)$,

con $latex 0 < x < 1$ y $latex u(0)= u(1)=0$. Particionamos el dominio en $latex n$ subintervalos introduciendo los puntos $latex x_i = ih$ con $latex h=frac{1}{n}$. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas (a partir de Taylor alrededor de $latex x_i$):

$latex frac{partial^2}{partial x^2} u(x_i) = frac{u_{i-1} -2 u_i + u_{i+1}}{h^2} – frac{h^2}{12} frac{partial^4}{partial x^4}u(xi_i)$

con $latex xi_i in (x_{i-1}, x_{i+1})$ por lo que:

$latex frac{v_{i-1}-2v_i+v_{i+1}}{h^2} = f_i$.

con un error de truncamiento local del orden de $latex O(h^2)$ y donde utilizamos la $latex v$ para hablar de la aproximación a la $latex u$. Si lo escribimos de manera matricial, tenemos:

$latex frac{1}{h^2} begin{bmatrix} -2 & 1 & & & & & \ 1 & -2 & 1 & & & & \ & 1 & -2 & 1 & & & \ & & cdot & cdot & cdot & & \ & & & cdot & cdot & cdot & \ & & & & 1 & -2 & 1 \ & & & & & 1 & -2end{bmatrix} begin{bmatrix} v_1 \ v_2 \ v_3 \ cdot \ cdot \ v_{n-2} \ v_{n-1} end{bmatrix} = begin{bmatrix} f_1 \ f_2 \ f_3 \ cdot \ cdot \ f_{n-2} \ f_{n-1} end{bmatrix}$

Existen otros esquemas de mayor orden para la aproximación, por ejemplo el esquema central de cuarto orden:

$latex frac{partial^2}{partial x^2} u(x_i) = frac{-u_{i-2}+16u_{i-1}-30u_{i}+16u_{i+1}-u_{i+2}}{12 h^2} + frac{h^4}{90} frac{partial^6}{partial x^6}u(xi_i)$

con $latex xi_i in (x_{i-2},x_{i+2})$.

En $latex 2D$ tenemos:

$latex u_{xx}+u_{yy} = f(x,y)$,

con $latex 0 < x < 1$, $latex 0 < y < 1$ y $latex u=0$ en el cuadrado unidad. Particionamos el dominio en $latex n_x times n_y$ subintervalos introduciendo los puntos $latex (x_i,y_j) = (ih_x, jh_y)$ con $latex h_x=frac{1}{n_x}$ y $latex h_y = frac{1}{n_y}$. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas:

$latex frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{h_x^2}+frac{v_{i,j-1}-2v_{i,j}+v_{i,j+1}}{h_y^2} = f_{i,j}$

con $latex v_{0,j} = v_{n_x,j} = v_{i,0} = u_{i,n_y} = 0$. Tenemos un error local de truncamiento del orden de $latex O(h_x^2+h_y^2)$.

Si hacemos $latex h_x = h_y = h$ nos queda una plantilla, que representa un operador que interactua localmente en la malla, como esta:

$latex A = frac{1}{h^2} begin{bmatrix} & 1 & \ 1 & -4 & 1 \ & 1 & end{bmatrix}$

En $latex 3D$ tenemos:

$latex u_{xx} + u_{yy} + u_{zz} = f(x,y,z)$,

con $latex 0 < x, y, z < 1$ y $latex u=0$ en la frontera del cubo unitario. Procediento de la misma manera que antes tenemos:

$latex frac{v_{i-1,j,k}-2v_{i,j,k}+v_{i+1,j,k}}{h_x^2}+frac{v_{i,j-1,k}-2v_{i,j,k}+v_{i,j+1,k}}{h_y^2} + frac{v_{i,j,k-1}-2v_{i,j,k}+v_{i,j,k+1}}{h_z^2}= f_{i,j,k}$

En este caso, la manera mas sencilla de representar el operador local es mediante el tensor:

$latex A = frac{1}{h^2} a_{ijk}$ con $latex a_{0jk}=begin{bmatrix} & & \ & 1 & \ & & end{bmatrix}, a_{1jk}=begin{bmatrix} & 1 & \ 1 & -6 & 1 \ & 1 & end{bmatrix}$ y $latex a_{2jk} = begin{bmatrix} & & \ & 1 & \ & & end{bmatrix}$

si $latex h_x = h_y = h_z = h$.

Tags: , ,

Aquí está el artículo donde aparece el nuevo esquema en el que el sistema se desacopla de manera jerárquica:

(1) Conocidas las cantidades hidrodinámicas conservadas, resolver:

$latex Delta X^i + frac{1}{3} mathcal{D}^i mathcal{D}_j X^j = 8 pi f^{ij} S_j^*$

para encontrar

$latex hat{A}^{ij} approx (LX)^{ij} = mathcal{D}^i X^j + mathcal{D}^j X^i – frac{2}{3} mathcal{D}_k X^k f^{ij}$.

(2) Resolver la ecuación:

$latex Delta psi = -2 pi psi^{-1} E^{*} – psi^{-7} frac{ f_{il} f_{jm} hat{A}^{lm} hat{A}^{ij} }{8}$

para encontrar $latex psi$, donde la unicidad local ahora esta garantizada. Podemos calcular $latex S^*$ de manear consistente.

(3) Resolver la ecuación:

$latex Delta(psi N) = 2 pi N psi^{-1} (E^* + 2 S^*) + N psi^{-7} frac{7 f_{il} f_{jm} hat{A}^{lm} hat{A}^{ij} }{8}$

para $latex N psi$, una ecuación lineal donde podemos aplicar el principio del máximo con lo que, con las codiciones de contorno apropiadas, se sigue la unicidad y existencia.

(4) Finalmente, resolver:

$latex Delta beta^i + frac{1}{3} mathcal{D}^i ( mathcal{D}_j beta^j ) = D_j( 2 N psi^{-6} hat{A}^{ij} )$

para encontrar $latex beta^i$.

Además, en este otro artículo, presentan una manera de reducir una ecuación elíptica vectorial, un complicado sistema acoplado de PDEs, a un conjunto de ecuaciones Poisson escalares desacopladas. Para el caso del shift, la $latex beta$ anterior, por ejemplo, en coordenadas esféricas, tendríamos:

(1) Resolver ecuación:

$latex Delta mu = mu_S$

que corresponde a la parte toroidal, para la resolución de la parte angular se introducen un potencial poloidal $latex eta$ y  un potencial toroidal $latex mu$ de manera que $latex boldsymbol{beta} = $, y está desacoplada del resto para obtener $latex mu$.

(2) Resolver la también desacoplada ecuación para la divergencia (de $latex boldsymbol{beta}$ respecto de la conexión plana $latex mathcal{D}$):

$latex Delta Theta = frac{3}{4} mathcal{D}_{hat{k}} S(boldsymbol{beta}^{hat{k}})$.

(3) Obtener $latex beta^r$ a partir de una de las siguiente ecuaciones:

(i) $latex frac{partial^2 beta^r}{partial r^2} + frac{4}{r} frac{partial beta^r}{partial r} + frac{2 beta^r}{r^2} + frac{1}{r^2}Delta_{theta varphi} beta^r = S(boldsymbol{beta})^r – frac{1}{3} frac{partial Theta}{partial r} + frac{2}{r} Theta$

(ii) $latex Delta chi = r S(boldsymbol{beta})^r – frac{r}{3} frac{partial Theta}{partial r} + 2 Theta $, donde $latex chi = r beta^r$

(4) Deducir $latex eta$ de una de las siguientes ecuaciones:

(i) $latex Delta_{theta varphi} eta = r Theta – r frac{partial beta^r}{delta r} – 2 beta^r$, que tiene la ventaja de que solo requiere una división por $latex -l (l+1)$ de los coeficientes de la expansión por armónicos esféricos pero la desventaja de que utiliza la derivada radial de $latex beta^r$ que puede tener problemas con el orden.

(ii) $latex Delta eta = eta_S – frac{2 beta^r}{r^2} – frac{1}{3} frac{Theta}{r}$, que requiere la resolución de otra ecuación de Poisson adicional.

Tags: , , , , ,