Ecuación de Poisson con condiciones de Dirichlet

Resumen: Se resuelve la ecuación de Poisson en un círculo con condición de contorno de tipo Dirichlet homogénea sobre toda la frontera.

Problema considerado

$\Omega\subset \mathbb{R}^2$ $f:\Omega \to \mathbb{R}$ $L^2(\Omega)$ $u:\Omega \to \mathbb{R}$ tal que

\begin{matrix} (P) & \begin{matrix} {\begin{cases} - Δ u = f & en Ω, \\ u = 0 & sobre \partial Ω . \end{cases} \end{matrix} \end{matrix}

Formulación variacional

Observación: todos los cálculos que se hacen en esta sección son puramente formales, es decir, se supone siempre la regularidad suficiente para que todos ellos estén justificados.

$\eqref{eq.poisson1}$ es el espacio de Sóbolev

\begin{matrix} H_{0}^{1} (Ω) = {v \in H^{1} (Ω) / v = 0 sobre \partial Ω} . \end{matrix}

$\eqref{eq.poisson1}$ $v\in H_0^1(\Omega)$ $\Omega$ :

\begin{matrix} \int_{Ω} (- Δ u) v d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω) . \end{matrix}

Integrando ahora por partes en el primer miembro, obtenemos

\begin{matrix} \int_{Ω} \nabla u \cdot \nabla v d x d y - \int_{\partial Ω} \frac{\partial u}{\partial n} v d σ = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω), \end{matrix}

$v$ $\partial\Omega$ , se obtiene:

\begin{matrix} \int_{Ω} \nabla u \cdot \nabla v d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω) . \end{matrix}

$\eqref{eq.poisson1}$ es, por tanto

\begin{matrix} (PV) & {\begin{cases} Hallar u \in H_{0}^{1} (Ω) tal que \\ \int_{Ω} \nabla u \cdot \nabla v d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω), \end{cases} \end{matrix}

$\nabla u \cdot \nabla v$ ,

\begin{matrix} {\begin{cases} Hallar u \in H_{0}^{1} (Ω) tal que \\ \int_{Ω} (\partial_{x} u \partial_{x} v + \partial_{y} u \partial_{y} v) d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω) . \end{cases} \end{matrix}

Problema aproximado

$P_1$ $\partial\Omega$ $\overline{\Omega}$ $\mathcal{T}_h = \{T_r\}^R_{r=1}$ $\overline{\Omega} = \bigcup_{r=1}^R T_r$ $\mathring T_r\neq \varnothing$ $r = 1,\dots, R$ $r\neq s$

\begin{matrix} T_{r} \cap T_{s} = {\begin{cases} un lado com \overset{´}{u} n a ambos tri \overset{´}{a} ngulos o \\ un v \overset{´}{e} rtice com \overset{´}{u} n a ambos tri \overset{´}{a} ngulos o \\ \emptyset \end{cases} \end{matrix}

Consideramos

\begin{matrix} W_{h} = {v_{h} | v_{h} \in C^{0} (\overset{―}{Ω}), v_{h} |_{T} \in P_{1} (T) \forall T \in T_{h}} \subset H^{1} (Ω), \end{matrix}

$W_h$ $H^1(\Omega)$ $\text{dim}\, W_h = M$ $M$ $\mathcal{T}_h$ . Consideramos también

\begin{matrix} V_{h} = {v_{h} | v_{h} \in W_{h}, v_{h} |_{\partial Ω} = 0} \subset H_{0}^{1} (Ω), \end{matrix}

$V_h$ $H_0^1(\Omega)$ $\text{dim}\, V_h = M_0$ $M_0$ $\mathcal{T}_h$ no $\partial\Omega$ $\eqref{eq.poisson1.var}$ es

\begin{matrix} (PV_{h}) & {\begin{cases} Hallar u_{h} \in V_{h} tal que \\ \int_{Ω} \nabla u_{h} \cdot \nabla v_{h} d x d y = \int_{Ω} f v_{h} d x d y \forall v_{h} \in V_{h} . \end{cases} \end{matrix}

$\{a_i\}_{i\in I}$ $\mathcal{T}_h$ $\{a_i\}_{i\in I_0}$ $I_0\subset I$ $\mathcal{T}_h$ $\partial\Omega$ $\text{card}(I) = M$ $\text{card}(I_0) = M_0$ .

$v_h\in W_h$ $\{v_h(a_i)\}_{i\in I}$ $v_h\in V_h$ $\{v_h(a_i)\}_{i\in I_0}$ .

$W_h$ $\{\varphi_i\}_{i\in I}$ definidas por

\begin{matrix} φ_{i} \in W_{h} \forall i = 1, \dots, M, φ_{i} (a_{j}) = δ_{i j} \forall i, j = 1, \dots, M . \end{matrix}

$\varphi_i$ $a_i$ $\mathcal{T}_h$ y solo es distinta de cero en los triángulos que contienen a dicho vértice.

$V_h$ $W_h$ $\{\varphi_i\}_{i\in I_0}$ . Luego,

\begin{matrix} v_{h} \in W_{h} ⟺ v_{h} = \sum_{i \in I} v_{h} (a_{i}) φ_{i}, \end{matrix}

\begin{matrix} v_{h} \in V_{h} ⟺ v_{h} = \sum_{i \in I_{0}} v_{h} (a_{i}) φ_{i} . \end{matrix}

$\{\varphi_i\}_{i\in I_0}$ $V_h$ el problema anterior se puede escribir

\begin{matrix} (1) & {\begin{cases} Hallar u_{h} \in V_{h} tal que \\ \int_{Ω} \nabla u_{h} \cdot \nabla φ_{i} d x d y = \int_{Ω} f φ_{i} d x d y \forall i \in I_{0} . \end{cases} \end{matrix}

$u_h$ $V_h$ :

\begin{matrix} u_{h} \in V_{h} ⟺ u_{h} = \sum_{j \in I_{0}} u_{h} (a_{j}) φ_{j}, \end{matrix}

de donde para $u_h$ $u_h(a_j)$ $j\in I_0$ .

$u_h$ $\eqref{eq.poisson1.var2}$ , el problema que tenemos que resolver se escribe:

\begin{matrix} (2) & {\begin{cases} Hallar u_{h} (a_{j})_{j \in I_{0}} tales que \\ \sum_{j \in I_{0}} u_{h} (a_{j}) \int_{Ω} \nabla φ_{j} \cdot \nabla φ_{i} d x d y = \int_{Ω} f φ_{i} d x d y \forall i \in I_{0} . \end{cases} \end{matrix}

Utilizando la siguiente notación

\begin{matrix} \begin{matrix} {\begin{cases} u_{j} = u_{h} (a_{j}), \\ a_{i j} = \int_{Ω} \nabla φ_{j} \cdot \nabla φ_{i} d x d y, \\ s_{i} = \int_{Ω} f φ_{i} d x d y, \end{cases} \end{matrix} \end{matrix}

$\eqref{eq.poisson1sl}$ se puede reescribir de la forma más compacta

\begin{matrix} (3) & {\begin{cases} Hallar (u_{j})_{j \in I_{0}} tales que \\ \sum_{j \in I_{0}} a_{i j} u_{j} = s_{i} \forall i \in I_{0}, \end{cases} \end{matrix}

$\eqref{eq.poisson1sl}$ sistema lineal $M_0$ ecuaciones.

Bloqueo de las condiciones de contorno

$(a_{ij})_{i,j\in I_0}$ $(a_{ij})_{i,j\in I}$ .

$~\eqref{eq.poisson1sl2}$ .

$\eqref{eq.poisson1sl2}$ es equivalente a

\begin{matrix} (4) & {\begin{cases} Hallar (u_{j})_{j \in I} tales que \\ {\begin{cases} \sum_{j \in I_{0}} a_{i j} u_{j} = s_{i} & i \in I_{0} \\ u_{i} = 0 & i \in I ∖ I_{0} . \end{cases} \end{cases} \end{matrix}

$u_h$ $\partial\Omega$ $M$ ecuaciones.

Este sistema, a su vez, es equivalente a

\begin{matrix} (5) & \sum_{j \in I} {\tilde{a}}_{i j} u_{j} = {\tilde{s}}_{i}, i \in I, \end{matrix}

donde

\begin{matrix} (6) & \begin{matrix} {\tilde{a}}_{i j} = {\begin{cases} V & si i = j \in I ∖ I_{0}, \\ a_{i j} & si no, \end{cases} {\tilde{s}}_{i} = {\begin{cases} V \cdot 0 \equiv 0 & si i = j \in I ∖ I_{0}, \\ s_{i} & si no, \end{cases} \end{matrix} \end{matrix}

$V$ muy grande $\approx 10^{30}$ ).

$i\in I\setminus I_0$ $i$ $\eqref{eq.poisson1sl3}$ :

\begin{matrix} (7) & \sum_{j \in I} {\tilde{a}}_{i j} u_{j} = \sum_{j \in I, j \neq i} a_{i j} u_{j} + {\tilde{a}}_{i i} u_{i} = {\tilde{s}}_{i} \Leftrightarrow \sum_{j \in I, j \neq i} a_{i j} u_{j} + V u_{i} = 0 \end{matrix}

$V$ $u_i\approx 0$ $i\in I\setminus I_0$ $u_i\approx 0$ $i\in I\setminus I_0$ $\eqref{eq.poisson1sl3}$ $i\in I_0$ $\eqref{eq.poisson1sl2}$ .

$\eqref{eq.poisson1sl3}$ $\eqref{sis:atilde}$ es el sistema que, realmente, se resuelve. Inicialmente, se construye la matriz y el segundo miembro del sistema "sin penalizar"

\begin{matrix} \sum_{j \in I} a_{i j} u_{j} = s_{i}, i \in I, \end{matrix}

$\eqref{sis:atilde}$ . El cálculo de la matriz y del segundo miembro de este sistema lineal se puede consultar en las notas sobre implementación del Método de los Elementos Finitos, aunque para resolver el sistema con FreeFEM no es necesario entrar en ese detalle.

Resolución con FreeFEM

$\Omega$ $\partial\Omega$ , viene definida por la curva

\begin{matrix} (8) & C = {(x, y) : x = \cos (t), y = \sin (t), 0 \leq t \leq 2 π} . \end{matrix}

$f(x,y) = xy$ .

Definición de la frontera del dominio

Las siguientes instrucciones definen la frontera del dominio, utilizando sus ecuaciones paramétricas y la dibujan. Obsérvese, que hemos asignado la etiqueta 1 con el comando label = 1 a toda la curva.


// Definicion de la frontera del dominio
border C(t = 0, 2*pi){x = cos(t); y = sin(t); label = 1;};
plot(C(50), wait = 1);

$C$ se han usado las variables reservadas x e y que representan, en FreeFEM, las dos variables espaciales.

Construcción del mallado

mesh $C$ discretizada con 60 segmentos sobre su contorno.


// Construccion de la triangulacion
mesh Th;
Th = buildmesh(C(60));
plot(Th, wait = 1);

Definición del espacio de elementos finitos

$W_h$ $\overline{\Omega}$ $T\in\mathcal{T}_h$ es un polinomio de grado menor o igual que 1.


// Definicion del espacio de EF
fespace Wh(Th, P1);

Definición del problema variacional

$f(x,y) = xy$ .


 // Definicion de la funcion del segundo miembro
func f = x*y;

ver aquí la técnica del bloqueo $u=0$ $\partial\Omega =C$ bloqueo $V_h$ $W_h$ .

En la definición del problema variacional con FreeFEM, la ecuación variacional se escribe en forma homogénea (segundo miembro = 0).

Además, los términos bilineales y lineales no deben estar bajo la misma integral. De hecho, para construir el sistema lineal, FreeFEM averigua qué integral contribuye a la forma bilineal, comprobando si están presentes ambos, la incógnita, que aquí se denota por u , y las funciones "test", denotadas aquí por v .


// Declaracion de las variables del problema (elementos de Wh)
Wh u, v;

// Definicion del problema variacional: 
problem MiPoisson(u, v) = 
        int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) 
       - int2d(Th)( f*v )
       + on(1, u = 0);

Aqui,

int2d es la función de FreeFEM que sirve para calcular integrales de superficie (en dominios 2D)
int2d(Th)( expresion )
es decir

\begin{matrix} \int_{Ω} expresion d x d y \end{matrix}

on es la función que impone la condición de Dirichlet (bloqueo):
```
on(ilabel, u = u0 )
```
bloquea ( a u0 ) el valor de la incógnita u en los nodos de la parte de frontera con etiqueta ilabel.
problem es un tipo de datos que permite definir la formulación variacional que se escribe en forma homogénea (segundo miembro = 0).
```
problem NombrePV(u, v) = ...
```
donde NombrePV es el nombre que se le da al problema variacional, en nuestro caso MiPoisson.

Resolución del problema y representación

Para resolver el problema, basta escribir el nombre del problema variacional:


// Resolucion
MiPoisson;

// Dibujo de la solucion
plot(u, wait = 1);

Complementos

$W_h$ usando el comando Wh.ndof .


// Definicion del espacio de EF
fespace Wh(Th, P1);
// Dimension de Wh
cout << "   " << endl;
cout << "Dimension del espacio de elementos finitos: " + Wh.ndof << endl;
cout << "   " << endl;

Las siguientes opciones de plot son interesantes:


// Opciones del dibujo
// plot( u, fill = true, nbiso = 40, value = true, dim = 3, wait = 1);

$u$ $W_h$ ${\cal T}_h$ . Este vector se denota u[]. Podemos asignar los valores de la solución a un vector.


// Recuperacion de los valores de la solucion en un vector
real[int] vv = u[];

cout << vv << endl;
// se obtiene lo mismo con
cout << u[] << endl;

Código completo del programa


// Definicion de la frontera del dominio
border C(t = 0, 2*pi){x = cos(t); y = sin(t); label = 1;};

// Construccion de la triangulacion
mesh Th;
Th = buildmesh(C(60));

// Definicion del espacio de elementos
fespace Wh(Th, P1);

// Definicion de la funcion del segundo miembro
func f = x*y;

// Declaracion de las variables del problema
Wh u, v;

// Definicion del problema variacional
problem MiPoisson(u, v) = 
        int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) 
    - int2d(Th)( f*v )
    + on(1, u = 0);

// Resolucion     
MiPoisson;

// Dibujo de la solucion
plot(u, wait = 1);

Ejercicios

Se proponen los siguientes ejercicios como variantes del problema:

$P_2$ $P_1$ , usando una triangulación grosera.
Usar mallados mas y menos finos.

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla