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 : \mathbb R^n \to \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 del problema 1D

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

1.2 Discretización por diferencias finitas

Sea una EDP definida en el dominio continuo \(x\in\Omega=[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} \} \]

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. \]

1.4 El operador discreto y sus propiedades

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

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

El apartado (a) de este lema implica que la matriz relacionada con \(\mathcal L_h\) es simétrica, y el apartado (b) implica que la matriz es definida positiva, luegos sus autovalores son positivos, existiendo solución y siendo única. Luego, el apartado (c) es la condición de estabilidad de la solución, y para problemas elípticos no depende de la discretización. El apartado (d) nos da una cota para \(\mathcal L_h\) y el apartado (e) nos da el espectro del operador.

Veamos una propiedad interesante que usaremos en varias demostraciones.

Ahora fijémonos en que si \(u_h, f_h\in V_h^0\) y definimos el operador \(\mathcal L_h\) tal que \(\mathcal L_h u_h = f_h\), entonces \(u_h = \mathcal L_h^{-1} f_h\), y aplicando el Lema 4c llegamos a: \[ ||\mathcal L_h^{-1} f_h||_{h,\infty} \le \frac18 ||f_h||_{h,\infty} \] Entonces, para puntos donde \(f_k \neq 0\), tenemos que \[ \frac{||\mathcal L_h^{-1} f_h||_{h,\infty}}{||f_h||_{h,\infty}} \le \frac18 \]

Este hecho nos inspira a definir el siguiente concepto:

Según esta definición, para el caso \(\mathcal L_h = -\Delta_h\), tendríamos \[ ||\mathcal L_h^{-1}||_{h, \infty} \le \frac18 \]

En la práctica queremos resolver \(u_h = \mathcal L^{-1} f_h\). El hecho de que el operador \(L_h^{-1}\) no mande funciones al infinito, i.e., mantiene a la función de entrada finita y controlada, nos garantiza la estabilidad del método.

Introduzcamos para el operador \(\mathcal L_h\) y el producto interno que hemos estado usando, la siguiente forma bilineal: \[ a_h (v_h, w_h) = (\mathcal L_h v_h, w_h)_h \]

Con este concepto vamos a definir la coercitividad:

La importancia de la coercitividad es que nos regala el siguiente resultado:

1.5 Análisis del error

Apreciemos que es posible relacionar de forma sencilla \(\tau_h\) con \(e_h\) de la siguiente forma: \[ \mathcal L_h e_h = L_h (u-u_h) = L_h u - L_h u_h = L_h u - f_h = \tau_h \] quedando que: \[ \boxed{ \mathcal L_h e_h = \tau_h} \]

1.6 Implementación computacional y verificación empírica

Para ilustrar los resultados teóricos obtenidos, en particular la resolución del sistema tridiagonal y el orden de convergencia \(p=2\), podemos implementar el operador \(-\Delta_h\) numéricamente.

El siguiente código en Python (utilizando numpy) construye la matriz del sistema \(\mathcal L_h u_h = f_h\) y resuelve el problema elíptico en 1D para un caso donde conocemos la solución analítica exacta: \(u(x) = \sin(2\pi x)\). Esta función satisface las condiciones de Dirichlet homogéneas \(u(0)=u(1)=0\), y su término fuente asociado es \(f(x) = -u''(x) = 4\pi^2\sin(2\pi x)\).

import numpy as np
import matplotlib.pyplot as plt


def finite_difference_elliptic_homogeneous_1D(f, N):
    h = 1/N

    x = np.linspace(0, 1, N+1)

    A = 2*np.diag(np.ones(N-1)) - np.diag(np.ones(N-2), 1) - np.diag(np.ones(N-2), -1)
    b = f(x[1:-1])

    u_num = np.zeros(N+1)
    u_num[1:-1] = h**2 * np.linalg.solve(A, b)

    return x, u_num

def finite_difference_elliptic_1D(f, a, b, N):
    h = 1/N
    def fh(x):
        fx = f(x)
        fx[0] += a/h**2
        fx[-1] += b/h**2
        return fx
    x, u_num = finite_difference_elliptic_homogeneous_1D(fh, N)
    u_num[0] = a
    u_num[-1] = b
    return x, u_num

if __name__ == "__main__":
    u = lambda x: np.sin(2*np.pi*x)
    f = lambda x: 4*np.pi**2*np.sin(2*np.pi*x)

    Ns = 2**np.arange(2, 10)
    hs = 1/Ns
    errors = np.empty_like(Ns, dtype=float)
    for i, N in enumerate(Ns):
        x, u_num = finite_difference_elliptic_1D(f, 1, 1, N)
        errors[i] = np.max(np.abs(u_num-u(x)))


    fig, ax = plt.subplots()

    ax.loglog(hs, errors, "or")
    ax.plot(hs, errors, "-r")

    ax.set_title("Error de discretización en función de la malla")
    ax.set_xlabel("h")
    ax.set_ylabel(r"$||e_h||_{h, \infty}$")

    p, logC = np.polyfit(np.log(hs), np.log(errors), 1)
    C = np.exp(logC)

    ax.text(0.87, 0.10, f" Pendiente: {p:0.2f}", transform=ax.transAxes, verticalalignment='top')

    ax.text(0.87, 0.05, r"$||e_h||_{h, \infty} \leq$" + f"{C:0.2f}" + r"$h^2$",
            transform=ax.transAxes, verticalalignment='top')
    
    plt.show()

Al ejecutar este script obtenemos la solución numérica, pero más valioso aún es que evaluamos empíricamente la cota del erro de discretización \(||e_h||_{h, \infty} \le Ch^2\). En una gráfica \(log-log\), la relación \(\text{Error} = Ch^p\) se transforma en la recta \(\log(\text{Error}) = p\log h + \log C\). El script hace un ajuste lineal con el método polyfit de numpy. La gráfica que nos devuelve el script es la siguiente:

Como vemos que la pendiente ajustada es \(2.02\), el código constata computacionalmente la convergencia cuadrática que demostramos de forma analítica en el Teorema 3.

1.7 Generalización para condiciones de contorno de Dirichlet no homogéneas

Hemos visto en el Apartado 1.3 la forma matricial del problema elíptico con condiciones de contorno homogéneas, en el que hemos usado el hecho de que \(u_0 = u_N = 0\) para reducir la matriz, pero en el fondo el sistema es: \[ \frac1{h^2} \begin{pmatrix} 1 & 0 & 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 \\ 0 & 0 & -1 & 2 & \cdots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots &\ddots & \vdots & \vdots &\vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 2 & -1 & 0 & 0 \\ 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_ 0 \\ u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-3} \\ u_{N-2} \\ u_{N-1} \\ u_N \end{pmatrix} = \begin{pmatrix} 0 \\ f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{N-3} \\ f_{N-2} \\ f_{N-1} \\ 0 \end{pmatrix} \]

Sea ahora el problema de Dirichlet no homogéneo: Encontremos \(u:\mathbb R \to \mathbb R\) tal que: \[ \begin{cases} -u''(x) = f & x\in [0, 1] \\ u(0) = a \\ u(1) = b \end{cases} \]

Si discretizamos nos queda resolver el siguiente sistema: \[ \frac1{h^2} \begin{pmatrix} 1 & 0 & 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\\ 0 & 0 & -1 & 2 & \cdots & 0 & 0 & 0 & 0\\ \vdots & \vdots & \vdots &\vdots &\ddots & \vdots & \vdots & \vdots &\vdots \\ 0 & 0 & 0 & 0 & \cdots & 2 & -1 & 0 & 0 \\ 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_ 0 \\ u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-3} \\ u_{N-2} \\ u_{N-1} \\ u_N \end{pmatrix} = \begin{pmatrix} \frac{a}{h^2} \\ f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{N-3} \\ f_{N-2} \\ f_{N-1} \\ \frac{b}{h^2} \end{pmatrix} \]

Una implementación en python que reutiliza la implimentación anterior (Apartado 1.6) es la siguiente:

import diferencias_finitas_eliptica_homogenea_1D as df_h
import numpy as np
import matplotlib.pyplot as plt

def finite_difference_elliptic_1D(f, a, b, N):
    h = 1/N

    def fh(x):
        fx = f(x)
        fx[0] += a/h**2
        fx[-1] += b/h**2
        return fx

    x, u_num = df_h.finite_difference_elliptic_homogeneous_1D(fh, N)
    u_num[0] = a
    u_num[-1] = b

    return x, u_num

if __name__ == "__main__":
    f = lambda x: 4*np.pi**2*np.sin(2*np.pi*x)
    N = 100
    x_nh, u_nh = finite_difference_elliptic_1D(f, -1, 1, N)
    x_h, u_h = df_h.finite_difference_elliptic_homogeneous_1D(f, N)

    fig, (ax1, ax2) = plt.subplots(1, 2)

    ax1.plot(x_h, u_h)
    ax1.set_title("Solución al problema homogéneo")
    ax1.set_xlabel("x")
    ax1.set_ylabel("u")

    ax2.plot(x_nh, u_nh)
    ax2.set_title("Solución al problema no homogéneo")
    ax2.set_xlabel("x")
    ax2.set_ylabel("u")

    plt.show()

1.8 Condición de contorno de Neumann

Enunciado: Encontrar \(u:\mathbb R \to \mathbb R\) tal que: \[ \begin{cases} -u''(x) = f & x\in[0, 1] \\ u'(0) = a &u(1) = b \end{cases} \]

Para abordar este caso no hacemos diferencias finitas hacia adelante para la condición de Neumann, ya que esto haría de nuestro esquema de orden \(O(h^2)\) un esquema de orden \(O(h)\). Para solucionar esto usaremos diferencias finitas centradas con un ‘nodo fantasma’: \[ \frac{u_1-u_{-1}}{2h} = a \implies u_{-1} = u_1 -2ah \]

Ahora podemos evaluar el operador discreto en el punto \(x_0\), que es un esquema de diferencias finitas centradas \(O(h^2)\), y sustituir por el resultado anterior: \[ \frac1{h^2}(-u_{-1}+2u_0-u_1) = f_0 \implies \frac1{h^2}(u_0-u_1) = \frac12f_0 - \frac ah \]

import diferencias_finitas_eliptica_homogenea_1D as df_h
import numpy as np
import matplotlib.pyplot as plt

def finite_difference_elliptic_neumann_1D(f, a, b, N):
    h = 1/N

    x = np.linspace(0, 1, N+1)

    A = 2*np.diag(np.ones(N)) - np.diag(np.ones(N-1), 1) - np.diag(np.ones(N-1), -1)
    B = f(x[0:-1])

    # Condición Neumann
    A[0, 0] = 1
    B[0] = B[0]/2 - a/h

    # Condición Dirichlet
    B[-1] += b/(h**2)

    u_num = np.zeros(N+1)
    u_num[0:-1] = h**2 * np.linalg.solve(A, B)
    u_num[-1] = b

    return x, u_num

if __name__ == "__main__":
    f = lambda x: 4*np.pi**2*np.sin(2*np.pi*x)
    x, u = finite_difference_elliptic_neumann_1D(f, 0, 1, 100)

    fig, ax = plt.subplots()

    ax.plot(x, u)
    ax.set_title("Solución al problema de contorno Neumann")
    ax.set_xlabel("x")
    ax.set_ylabel("u")

    plt.show()

Es muy importante resaltar que siempre debe haber una condición de Dirichlet. Supongamos el problema elíptico unidimensional \(-u''(x)=f\) con las condiciones de Neumann \(u'(0) = a\) y \(u'(1) = b\), entonces la función \(v(x)=u(x)+C\) también es solución de la misma EDP. Es decir, perdemos la unicidad de la solución, y obtenemos una familia de soluciones paralelas. De hecho, si integramos la EDP: \[ -\int_0^1 u''(x)dx = \int_0^1 f(x)dx \implies -(u'(1)-u'(0)) = \int_0^1 f(x)dx \implies \boxed{a-b = \int_0^1 f(x)dx} \] Si las condiciones de contorno no cumplen la condición anterior con la función fuente, perdemos la existencia de la solución.

1.9 Caso bidimensional

Enunciado: Encontrar \(u: \mathbb R^2 \to \mathbb R\) tal que: \[ \begin{cases} -\nabla^2 u(x, y) = 1 &(x, y) \in [0, 1]\times [0, 1] \\ u(0, y) = u(1, y) = 0 \\ u(x, 0) = u(x, 1) = 0 \end{cases} \]

Discretizamos el espacio de forma que \(x_i = ih\) e \(y_j = jh\), de modo que: \[ u_{i,j} = u(x_i, y_i) = u(ih, jh) \]

Usaremos el resultado unidimensional para definir el operador bidimensional: \[ \begin{align} \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} &\approx \frac{u_{i-1, j} - 2u_{i,j} + u_{i+1, j}}{h^2} + \frac{u_{i, j-1} - 2u_{i,j} + u_{i, j+1}}{h^2} \\ &= \frac{u_{i-1, j} + u_{i+1, j} + u_{i, j-1} + u_{i, j+1} - 4u_{i,j}}{h^2} \end{align} \]

Para poder representar el sistema de ecuaciones asociado a la discretización, primero debemos idear una forma de ‘aplanar’ la matriz, i.e, definir una \(k\) tal que \(u_{i, j} = u_k\).

(diagrama)

Por observación podemos identificar que \(k=(j-1)(N-1)+i\), pero necesitamos la relación inversa. Esta la obtenemos observando que: \[ \begin{align} j = 1 \implies &k \in \{1,\:2,\:\cdots,\:N-1\} \\ j = 2 \implies &k \in \left\{(N-1)+1,\:(N-1)+2,\:\cdots,\:2(N-1)\right\}\\ j = 3 \implies &k \in \{2(N-1)+1,\:2(N-1)+2,\:\cdots,\:3(N-1)\} \end{align} \]

Donde podemos identificar que: \[ \frac{k}{N-1} \in \left(\: j-1,\:j\:\right] \implies \boxed{ j = \left\lceil \frac{k}{N-1} \right\rceil } \]

Y a partir de \(j\) podemos despejar \(i\): \[ k=(j-1)(N-1)+i \implies \boxed{i = k - (j-1)(N-1)} \]

Fijémonos en que hemos ‘aplanado’ la solución discreta, pasando de una matriz \((N-1)\times (N-1)\) a un vector columna de \((N-1)^2\) elementos. Por tanto, la matriz asociada al sistema lineal será una matriz \((N-1)^2\times(N-1)^2\). Para construir el sistema lineal que resuelva la discretización, tomaremos una filosofía de bloques. Fijémonos de qué forma se relacionan los elementos de una fila con los otros elementos de la misma fila: \[ A = \begin{pmatrix} 4 & -1 & 0 & \cdots & 0 & 0 & 0 \\ -1 & 4 & -1 & \cdots & 0 & 0 & 0 \\ 0 & -1 & 4 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 4 & -1 & 0 \\ 0 & 0 & 0 & \cdots & -1 & 4 & -1 \\ 0 & 0 & 0 & \cdots & 0 & -1 & 4 \end{pmatrix}_{(N-1)\times(N-1)} \]

Y ahora de qué forma se relacionan los elementos de una fila con los elementos de la fila contigua (siguiente u anterior): \[ B = \begin{pmatrix} -1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & -1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & -1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & -1 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & -1 \end{pmatrix}_{(N-1)\times(N-1)} \] Finalmente la matriz del sistema será la matriz compuesta por estos bloques: \[ M = \begin{pmatrix} A & B & 0 & \cdots & 0 & 0 & 0 \\ B & A & B & \cdots & 0 & 0 & 0 \\ 0 & B & A & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & A & B & 0 \\ 0 & 0 & 0 & \cdots & B & A & B \\ 0 & 0 & 0 & \cdots & 0 & B & A \end{pmatrix}_{(N-1)^2\times(N-1)^2} \]

Para visualizar por qué esto es así, debemos tener en cuenta que un elemento en la columna \(i\) de la fila \(j\) tendrá relación, fuera de su propia fila, con los elementos \((i, j-1)\) y \((i, j+1)\), que en la forma ‘aplanada’ se encuentran a una distancia \(N-1\) de nuestro elemento de referencia. Esto hace los elementos de la diagonal del bloque anterior y siguiente sean \(-1\), formando el patrón que hemos visto.

Por ejemplo, para el caso en el que \(N=4\): \[ M= \begin{pmatrix} 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 4 & 0 & 0 & -1 & 0 & 0 & 0 \\ -1 & 0 & 0 & 4 & -1 & 0 & -1 & 0 & 0 \\ 0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0 \\ 0 & 0 & -1 & 0 & -1 & 4 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 & 0 & 0 & 4 & -1 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 \\ \end{pmatrix} \]

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from matplotlib import cm


def finite_difference_elliptic_homogeneous_2D(f, N):
    h = 1 / N

    # Definimos las submatrices
    A = sp.diags([-1, 4, -1], [-1, 0, 1], shape=(N - 1, N - 1), dtype=np.float64)
    B = -sp.eye(N - 1)

    # Definimos la distribución de las submatrices
    F_A = sp.eye(N - 1)
    F_B = sp.diags([1, 1], [-1, 1], shape=(N - 1, N - 1), dtype=np.float64)

    # Usamos la delta de Kronecker para montar la matriz
    M = sp.kron(F_A, A) + sp.kron(F_B, B)

    # Comprimimos a formato CSR (para cálculos)
    M = M.tocsr()

    # Calculamos los términos independientes
    points = np.linspace(0, 1, N + 1)
    inner_points = points[1:-1]

    X, Y = np.meshgrid(inner_points, inner_points, indexing='ij')

    b = h ** 2 * f(X, Y).ravel()

    # Resolvemos y montamos la matriz solución
    u_vec = sp.linalg.spsolve(M, b)

    U_inner = u_vec.reshape((N - 1, N - 1))

    U = np.zeros([N + 1, N + 1])
    U[1:-1, 1:-1] = U_inner

    X, Y = np.meshgrid(points, points)

    return X, Y, U


if __name__ == "__main__":
    f = lambda x, y: 8 * np.pi ** 2 * np.sin(2 * np.pi * x) * np.sin(2 * np.pi * y)
    X, Y, U = finite_difference_elliptic_homogeneous_2D(f, 100)

    fig = plt.figure(figsize=(10, 7))
    ax = fig.add_subplot(111, projection='3d')

    surf = ax.plot_surface(X, Y, U, cmap=cm.viridis,
                           linewidth=0, antialiased=True)

    ax.set_title(r'Solución de la Ecuación de Poisson $\nabla^2 u = f$')
    ax.set_xlabel('Eje X')
    ax.set_ylabel('Eje Y')
    ax.set_zlabel(r'$u(x,y)$')
    fig.colorbar(surf, shrink=0.5, aspect=5)
    plt.show()

2. Problemas parabólicos

El ejemplo por excelencia de un problema parabólico es la ecuación del calor con condiciones de contorno de Dirichlet homogéneas y una condición inicial:

Enunciado: encuentra \(u : \mathbb [0, T] \times R^n \to \mathbb R\) que satisfaga: \[ \left\{ \begin{align*} &\frac{\partial u}{\partial t} - \alpha \nabla^2 u = f(t, \vec x) && \forall t \in (0, T], \forall \vec x \in\Omega \subseteq \mathbb{R}^n\\ &u(\vec x, t) = 0 && \forall t\in (0, T], \forall \vec x \in \Gamma = \partial \Omega \\ &u(\vec x, 0) = u_0(\vec x) && \forall \vec x \in \Omega \end{align*} \right. \]

donde \(\alpha > 0\) es la difusividad térmica.

2.1 Solución analítica del problema 1D

Resolvamos el caso particular cuando \(\Omega = [0, 1]\) y no haya término fuente: \[ \left\{ \begin{align*} &\frac{\partial u}{\partial t}-\frac{\partial^2 u}{\partial x^2} = 0 && \forall t \in (0, T], \forall x \in [0, 1]\\ &u(0, t) = u(1, t) = 0 && \forall t\in (0, T] \\ &u(x, 0) = u_0(x) && \forall x \in [0,1] \end{align*} \right. \]

Usaremos el método de separación de variables, que consiste en suponer que la solución es de la siguiente forma: \[u (x, t) = X(x) T(t)\]

Por tanto, sustituyendo en nuestra ecuación: \[ X(x)\frac{dT}{dt} - T(t)\frac{d^2 X}{dx^2} = 0 \]

Reordenamos para que cada término dependa de una única variable: \[ \frac1T\frac{dT}{dt} = \frac1X\frac{d^2 X}{dx^2} = -\lambda \] que hemos igualado a una constante debido a que cada término depende de una única variable, y la igualdad se cumple para todo el dominio, lo que solo puede pasar si cada término es constante. Resolviendo estas dos EDOs para las condiciones de contorno dadas, nos quedan las siguientes infinitas soluciones para cada \(k\in\mathbb Z\): \[ \left\{ \begin{align*} &X_k (x) = \sin(k\pi x) \\ &T_k (t) = e^{-\left(k\pi\right)^2 t} \end{align*} \right. \]

Quedándonos las siguientes infinitas soluciones: \[ u_k(x, t) = \sin\left(k\pi x\right) e^{-\left(k\pi\right)^2 t} \]

2.2 Discretización por diferencias finitas

Vamos a discretizar el dominio \([0, T]\times[0, 1]\) en \(N\) intervalos espaciales y \(M\) intervalos temporales, de forma que: \[ \left\{ \begin{align*} &x_i = x_0 + ih &i=0,1,\cdots, N\\ &t^k = t_0 + k \Delta t &k=0,1,\cdots, M \end{align*} \right. \]

Lo haremos en dos pasos.

• Discretizamos el espacio y aplicamos diferencias finitas. \[ \frac{\partial u_i}{\partial t} = \frac{u_{i-1}-2u_i+u_{i+1}}{h^2} \]
• Aproximamos el anterior sistema de EDOs usando Euler explícito: \[ u_i^{k+1} = u_i^k + \Delta t \left( u_{i-1} -2u_i +u_{i+1} \right) \]

Si ahora despejamos el único término en \(k+1\), tenemos: \[ u_i^{k+1} = \lambda u_{i+1}^k + (1-2\lambda) u_i^k + \lambda u_{i-1}^k \] donde hemos definido \(\lambda = \frac{\Delta t}{h^2}\).