Métodos Numéricos para EDPs

David Vereb

Resumen

Estos apuntes son una síntesis personal de las clases impartidas por Albert Oliver Serra. Cualquier error u omisión son mi responsabilidad y no reflejan necesariamente el contenido oficial del curso.

Parte I: Métodos de Diferencias Finitas

1. Problemas elípticos

El ejemplo por excelencia de un problema elíptico es la Ecuación de Poisson con condiciones de contorno de Dirichlet homogéneas:

Enunciado: encuentra \(u(\vec x)\in\mathbb R\) que satisfaga: \[ \begin{cases} -\nabla^2 u = f(\vec x) & \forall \vec x \in\Omega \subseteq \mathbb{R}^n \\ u(\vec x) = 0 & \forall \vec x \in \Gamma = \partial \Omega \\ \end{cases} \]

1.1 Solución analítica

Resolvamos el caso particular cuando \(\Omega = [0, 1]\), en el que \(u\in C^2([0, 1])\), y: \[ \begin{cases} -u'' = f(x) & \forall x\in[0, 1] \\ u(0) = u(1) = 0 \end{cases} \]

Si integramos y definimos \(F(x)\) tal que \(F'(x) = f(x)\): \[ u' = -F(x) + C_1 \]

Como F(x) no es única, forzamos que \(F(x) = \int_0^x f(s)ds\), luego \(F(0) = 0\). Si volvemos a integrar: \[ u(x) = -\left(\int_0^x F(s)\,ds - C_1 x + C_2 \right)\] con \(C_1, C_2\in\mathbb{R}\).

Aplicando las condiciones de contorno: \[ \begin{cases} C_1 = \int_0^1 F(s)\,ds \\ C_2 = 0 \end{cases} \]

Reuniendo estas condiciones nos queda: \[u(x) = -\int_0^x F(s)\,ds + x\int_0^1 F(s)\,ds\]

Considerando la definición de \(F(s)\), podemos reescribir la expresión anterior: \[u(x) = \int_0^x s(1-x)f(s)\,ds + \int_x^1 x(1-s) f(s)\,ds\]

Y si definimos la función de Green: \[ G(x, s) = \min(x, s)[1-\max(x,s)] = \begin{cases} s(1-x) & \text{si } 0\le s\le x \\ x(1-s) & \text{si } x\le s\le 1 \end{cases} \]

Podemos finalmente escribir la solución como: \[ u(x) = \int_0^1 G(x, s) f(s)\,ds \]

Propiedades de la solución:

Veamos la demostración de esta última propiedad: \[ |u| \le \int_0^1 G(x,s) |f(s)|\,ds \le ||f||_\infty \int_0^1 G(x,s)\,ds = \frac12 x(1-x) ||f||_\infty \]

y como el máximo de \(x(1-x)\) está en \(x=\frac12 \Rightarrow x(1-x)=\frac14\), se demuestra lo propuesto.

1.2 Discretización del dominio de una EDP

Sea una EDP definida en el dominio continuo \(x\in\Omega\subseteq[0,1]\). Definimos el dominio discreto \(\Omega_h\) de \(N+1\) puntos y \(N\) intervalos, caracterizado por el parámetro \(h=\frac1N\) que representa la longitud de un intervalo.

Así, los puntos de la discretización son los siguientes: \[ \Omega_h = \{x_i=ih\::\:i=0,1,\dots,N\in\mathbb{N} \} \]

Definición: En \(\Omega_h\) tenemos que aproximar \(u''\) con diferencias finitas, por ejemplo usando un esquema centrado: \[ u''(x)\approx \frac{u(x-h)-2u(x)+u(x+h)}{h^2} := \Delta_h u(x) \] donde \(\Delta_h\) es el operador de segunda diferencia finita.

Lema: Sea \(v\in C^4([0, 1])\), entonces \(\forall x\in[h, 1-h]\), se tiene que: \[ \boxed{ |v''-\Delta_h v|| \le \frac{h^2}{12}||v^{(iv)}||_\infty } \]

Demostración del lema: aplicando el Teorema de Taylor a \(v(x+h)\) y \(v(x-h)\) obtenemos: \[ \begin{align*} v(x+h) = v(x) + hv'(x) + \frac1{2!} h^2 v''(x) + \frac1{3!} h^3 v'''(x) + \frac1{4!} h^4 v^{(iv)}(\xi_{+}) \\ v(x-h) = v(x) - hv'(x) + \frac1{2!} h^2 v''(x) - \frac1{3!} h^3 v'''(x) + \frac1{4!} h^4 v^{(iv)}(\xi_{-}) \end{align*} \] donde \(\xi_{+}\in[x, x+h]\) y \(\xi_{-}\in[x-h, x]\).

Si sumamos ambas expresiones, nos queda: \[ v(x-h)+v(x+h) = 2v(x)+h^2v''(x)+\frac1{4!}\left( v^{(iv)}(\xi_+) + v^{(iv)} (\xi_-)\right) \]

Podemos reorganizar la ecuación de forma que identifiquemos el operador de segunda diferencia finita: \[ \Delta_h v(x) = \frac{v(x-h)-2v(x)+v(x+h)}{h^2} = v''(x)+\frac{1}{4!}h^2\left( v^{(iv)}(\xi_+) + v^{(iv)}(\xi_-) \right) \]

Reorganizando términos, tomando valor absoluto y aplicando desigualdad triangular, llegamos a: \[ |v''(x) - \Delta_h v(x)| \le \frac1{4!}h^2\left( |v^{(iv)}(\xi_+)| + |v^{(iv)}(\xi_-)| \right) \]

Sea \(\xi\in\{\xi_+, \xi_-\}\) tal que \(v^{(iv)}(\xi) = \max\left(v^{(iv)}(\xi_+), v^{(iv)}(\xi_-)\right)\), entonces: \[ |v''(x) - \Delta_h v(x)| \le \frac{2}{4!}h^2|v^{(iv)}(\xi)|\le \frac{h^2}{12}||v^{(iv)}(x)||_\infty \] Quedando demostrado nuestro lema.

1.3 Solución de la discretización

Sea la siguiente EDP con condiciones de contorno: \[ \begin{cases} -u''= f(x) & x\in[0, 1] \\ u(0) = u(1) = 0 \end{cases} \]

Si discretizamos el dominio a \(\Omega_h\) con \(N\) intervalos: \[ \begin{cases} -\Delta_h u(x) = f(x) \\ u(0)=u(1)=0 \end{cases} \]

Lo que se traduce en el siguiente sistema de ecuaciones: \[ \left\{ \begin{align*} &-\frac{ u_{i-1} - 2u_i + u_{i+1}}{h^2} = f_i & i=1,\dots, N-1 \\ &u_0 = u_N = 0 \end{align*} \right. \]

El problema se reduce a resolver este sistema, que en forma matricial es una matriz tridiagonal simétrica: \[ \frac1{h^2} \begin{pmatrix} 2 & -1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & \cdots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots &\vdots &\ddots & \vdots & \vdots & \vdots &\vdots \\ 0 & 0 & 0 & 0 & \cdots & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-3} \\ u_{N-2} \\ u_{N-1} \end{pmatrix} = \begin{pmatrix} f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{N-3} \\ f_{N-2} \\ f_{N-1} \end{pmatrix} \]

1.4 Análisis del error

Definición: Llamamos función grid a una función \(\omega_h\) cuyo dominio es el espacio discreto \(\Omega_h\), i.e., \(\omega_h: \Omega_h \to \mathbb R\). El conjunto de todas las funciones grid lo denotamos \(V_h\), y definimos el siguiente subespacio: \[ V_h^0 = \{ \omega_h\in V_h : \omega_h(0) = \omega_h(1) = 0 \} \]

Definición: Definimos el operador lineal discreto \(\mathcal{L}_h : V_h^0 \to V_h^0\) tal que, para todo \(\omega\in V_h^0\): \[ \left\{ \begin{align*} &\mathcal{L}_h[\omega_h](x_i) = -\Delta_h \omega_h(x_i) & i=1,\dots, N-1 \\ &\mathcal{L}_h[\omega_h](x_0) = \mathcal{L}_h[\omega_h](x_N) = 0 \end{align*} \right. \]

Por lo tanto, podemos definir el problema como encontrar \(u_h\in V_h^0\) tal que: \[\boxed{\mathcal L u_h = f_h}\] con \(f_h = \left. f \right|_{\Omega_h}\).

Definición: Definimos la función discreta de Green \(G^k\in V_h^0\), para \(k=1,\dots, N-1\), como la solución de: \[ \mathcal L_h G^k = \delta^k \hspace{1cm} k=1, \dots, N-1 \\ \] donde \(\delta^k\in V_h^0\) tal que: \[\delta^k(x_i) = \delta(k-i) = \begin{cases} 0 & \text{si } i\neq k \\ 1 & \text{si } i=k\end{cases}\]

Lema: Se cumple la siguiente relación entre las funciones de Green discreta y continua. \[ G^k(x_i) = h G(x_i, x_k)\] para \(i=1,\dots, N-1\).

Demostración del lema: Sea \(x_k \in \Omega_h\).

El primer paso será mostrar que \(\mathcal L_h G(\,\cdot\,, x_k)(x_i) = \frac1h \delta^k(x_i)\). Distingamos dos casos:

Hemos visto que \(\mathcal L_h G(\,\cdot\,, x_k)(x_i) = \frac1h \delta^k(x_i)\), luego: \[ \mathcal L_h (hG(\,\cdot\,, x_k)) = \delta^k = \mathcal L_h G^k \] Como la matriz asociada al operador es definida positiva, entonces el operador es invertible, y en consecuencia inyectiva. Esto prueba nuestro lema.

1.5 Propiedades algebraicas y estabilidad

Definición: Definimos el producto interno en el espacio \(V_h^0\) de dos vectores \(v_h, w_h \in V_h^0\) como: \[ (v_h, w_h)_h = h\sum_{i=1}^{N-1} v_i w_i \] donde \(v_i = v_h(x_i)\) y \(w_i = w_h(x_i)\).

Definición: Definimos la norma infinito de un vector \(v_h\in V_h^0\) como: \[ ||v_h||_{h,\infty} = \max_{0\le i\le N}|v_h(x_i)| \]

A continuación se presenta un lema de propiedades del operador \(\mathcal L_h\):

Lema. Se cumplen las siguiente propiedades:

  1. Simetría: El operador \(\mathcal L_h\) es simétrico para todo \(v_h, w_h\in V_h^0\), i.e.: \[ (\mathcal L_h v_h, w_h)_h = (v_h, \mathcal L_h w_h)_h \]

  2. Definido positivo: El operador \(\mathcal L_h\) es definidio positivo, i.e., para todo \(v_h\in V_h^0\), se tiene que: \[ \left\{ \begin{align*} &(\mathcal L_h v_h, v_h)_h \ge 0 \\ &(\mathcal L_h v_h, v_h)_h = 0 \Leftrightarrow v_h = 0 \end{align*} \right. \]

  3. Estabilidad: Si \(\mathcal L_h u_h = f_h\), entonces: \[\boxed{||u_h||_{h,\infty} \le \frac18 ||f_h||_{h,\infty}}\]

  4. Acotación: El operador \(\mathcal L_h\) está acotado, i.e., para todo \(v_h, w_h\in V_h^0\), se tiene que: \[ (\mathcal L_h v_h, w_h)_h \le 4 ||v_h||_h ||w_h||_h \]

  5. Espectro: Los autovalores de \(\mathcal L_h\) están definidos, para cada \(k=1,\cdots,N-1\), como: \[ \lambda_k = \frac{4}{h^2} \sin^2\left(\frac{kh\pi}2\right) \] Y sus respectivos autovectores son \(\vec z_h^k\), cuya componente \(j\)-ésima es: \[ z_j^k = \sqrt{2h} \sin(jkh\pi) \]

La definida positividad (b) es la clave para la existencia y unicidad de la solución, ya que implica que el operador \(\mathcal L_h\) es inyectivo y, al estar en un espacio de dimensión finita, biyectivo; esto asegura que para cualquier función de carga \(f_h\), existe una única solución \(u_h\). Por otro lado, la estabilidad (c) asegura que el problema está bien definido, pues garantiza que la solución depende de forma continua de los datos de entrada, impidiendo que pequeños errores de redondeo o perturbaciones en \(f\) se amplifiquen sin control.

  1. Demostración simetría:

La demostración se basa en ajustar los límites de las sumas para reordenarlas como nos interesa. Usaremos el hecho de que \(v_h,w_h\in V_h^0\), luego \(v_0=w_0=v_N=w_N=0\), y el desarrollo es el que sigue: \[ \begin{align*} (\mathcal L_h v_h, w_h)_h &= h\sum_{i=1}^{N-1} \frac1{h^2} (-v_{i-1}+2v_i-v_{i+1})w_i = \\ &= -\frac1h \sum_{i=1}^{N-1}[(v_{i+1}-v_i)-(v_i-v_{i-1})]w_i = \\ &= -\frac1h \left[\sum_{i=1}^{N-1} (v_{i+1}-v_i)w_i - \sum_{i=0}^{N-2}(v_{i+1}-v_i)w_{i+1}\right] = \\ &= -\frac1h \left[\sum_{i=1}^{N-1} (v_{i+1}-v_i)w_i - \sum_{i=1}^{N-1}(v_{i+1}-v_i)w_{i+1} +v_1w_1 \right]= \\ &= \frac1h \left[\sum_{i=1}^{N-1} (v_{i+1}-v_i)(w_{i+1}-w_i)+v_1 w_1\right] = \\ &= \frac1h \left[\sum_{i=1}^{N-1} v_{i+1}(w_{i+1}-w_i) - \sum_{i=1}^{N-1}v_i(w_{i+1}-w_i) + v_1 w_1 \right] = \\ &= \frac1h \left[\sum_{i=2}^N v_i(w_i-w_{i-1}) - \sum_{i=1}^{N-1}v_i (w_{i+1}-w_i)+v_1 w_1 \right] = \\ &= \frac1h \left[-\sum_{i=1}^{N-1} v_i[(w_{i+1}-w_i)-(w_i-w_{i-1})] - v_1w_1 + v_1w_1 \right] = \\ &= -\frac1h \sum_{i=1}^{N-1} v_i[(w_{i+1}-w_i)-(w_i-w_{i-1}) = \\ &= h\sum_{i=1}^{N-1} v_i \frac1{h^2}(w_{i+1}-2w_i+w_{i-1}) = (v_h, \mathcal L_h w_h)_h \end{align*} \] Que es lo que queríamos demostrar.

  1. Demostración definido positivo:

Del desarrollo de la demostración de (a), para todo \(v_h\in V_h^0\), se tiene que: \[ (\mathcal L_h v_h, v_h)_h = \frac1h \left[\sum_{i=1}^{N-1} (v_{i+1}-v_i)^2 + v_i^2 \right] \ge 0 \] Además, si \((\mathcal L_h v_h, v_h) = 0\), entonces \(v_{i+1}=v_i\), y como \(v_1=0\), por inducción \(v_i=0\) para \(i=1,\cdots,N-1\). El recíproco también es cierto, ya que \(v_i=0\) para \(i=1,\cdots,N-1\) implica que \((\mathcal L_h v_h, v_h) = 0\).

  1. Demostración estabilidad: Hemos demostrado que el operador \(\mathcal L_h\) es definido positivo, luego es invertible. Por tanto, se cumple que: \[ \mathcal L_h v_h = f_h \Rightarrow u_h = \mathcal L_h^{-1} f_h \]