Mauricio Villa Amaya, Alejandro Sánchez Yali

Introducción

Considere el problema de valor inicial (PVI) dado por la ecuación de difusión en dos dimensiones: \begin{equation}\label{eq: difusión} \partial_{t}u(x,\,y,\,t)=\alpha \nabla^2 u(x,\,y,\,t), \end{equation} con condición inicial sobre el rectángulo $\Omega=[a,\, b]\times [c,\, d]$ dada por $u(x,y, 0)=f(x,\, y)$ y la condición de frontera de Dirichlet $u(x,\, y,\, t)=g(x,\,y,\,t)$ con $(x,\,y)\in \partial \Omega$ para todo $t>0$. El objetivo es solucionar la ecuación de difusión por los métodos: Forward-time central-space (FTCS) y backward-time central-space (BTCS).

Sean $I_{_l}= \{0,\, \dots,\, l\}$ y $I_{_m}=\{0,\,\dots,\, m\}$, y la partición, sobre el conjunto $\Omega$, $P_{_{l,m}}=\{(x_{_j},\, y_{_k}): (j,\, k)\in I_{_l}\times I_{_m}\}$, donde $x_{_j}= a+ j\Delta x$ con $\Delta x= \frac{b-a}{l}$ e $y_{_k}=c+k\Delta y$ con $\Delta y = \frac{d-c}{m}$. Observe que los números $l$ y $m$ se interpretan como el número de subintervalos en los que se dividen $[a,\, b]$ y $[c,\, d]$ respectivamente. Para la variable $t$, considere los $t_{_n}=n\Delta t$. Además para los puntos de la partición se empleará la notación $u^{n}_{_{j,k}} = u(a+j\Delta x,\,b+k\Delta y,\, n\Delta t)$.

Se definen sobre $u^{n}_{_{j,k}}$ las funciones \begin{eqnarray}\label{eq:suc-ant} s_{_j}(u_{_{j,k}}^n)=u_{_{j+1,k}}^n, \quad s_{_k}(u_{_{j,k}}^n)=u_{_{j,k+1}}^n, \quad s_{_n}(u_{_{j,k}}^n)=u_{_{j,k}}^{n+1},\\ a_{_j}(u_{_{j,k}}^n)=u_{_{j-1,k}}^n, \quad a_{_k}(u_{_{j,k}}^n)=u_{_{j,k-1}}^n, \quad a_{_n}(u_{_{j,k}}^n)=u_{_{j,k}}^{n-1}, \end{eqnarray} que satisfacen las siguientes propiedades:

1. $s_{\circ}(bu_{_{j,k}}^n+cu_{_{r,s}}^t)=bs_{\circ}(u_{_{j,k}}^n)+cs_{\circ}(u_{_{r,s}}^t)$,
2. $a_{\circ}(bu_{_{j,k}}^n+cu_{_{r,s}}^t)=ba_{\circ}(u_{_{j,k}}^n)+ca_{\circ}(u_{_{r,s}}^t)$,
3. $(s_{\circ}+s_{\bullet})(u_{_{j,k}}^n)=s_{\circ}(u_{_{j,k}}^n)+s_{\bullet}(u_{_{j,k}}^n)$,
4. $(s_{\circ}+a_{\bullet})(u_{_{j,k}}^n)=s_{\circ}(u_{_{j,k}}^n)+a_{\bullet}(u_{_{j,k}}^n)$,
5. $(a_{\circ}+a_{\bullet})(u_{_{j,k}}^n)=a_{\circ}(u_{_{j,k}}^n)+a_{\bullet}(u_{_{j,k}}^n)$,

donde $b, c$ son números reales. Las funciones $s_{\circ}$ y $a_{\circ}$ las denominaremos sucesor y antecesor de $\circ$ respectivamente.

Método explícito FTCS

Al aplicar el método FTCS en la ecuación \eqref{eq: difusión} se obtiene:
\begin{equation}\label{eq:ftcs} \dfrac{u_{_{j,k}}^{n+1} - u_{_{j,k}}^{n}}{\Delta t}\approx \alpha \left( \dfrac{u_{_{j+1,k}}^{n}-2u_{_{j,k}}^{n}+u_{_{j-1,k}}^{n}}{\Delta x^{2}} + \dfrac{u_{_{j,k+1}}^{n}-2u_{_{j,k}}^{n}+u_{_{j,k-1}}^{n}}{\Delta y^{2}} \right). \end{equation}
Empleando las funciones sucesor y antecesor en la ecuación \eqref{eq:ftcs}, el PVI discretizado toma la forma:

\begin{equation}\label{eq:ftcs2} \left\{\begin{array}{cc} \dfrac{(s_{_n}-\mathrm{id})}{\Delta t}(u_{_{j,k}}^{n})\approx \alpha\left(\dfrac{s_{_j}-2\mathrm{id}+a_{_j}}{\Delta x^{2}} + \dfrac{s_{_k}-2\mathrm{id}+a_{_k}}{\Delta y^{2}}\right)(u_{_{j,k}}^{n}) & \qquad \forall (x_{_j},\, y_{_k})\in P(\mathring{\Omega}),\\ u_{_{j,k}}^{0}=f_{_{j,k}} & \qquad \forall (x_{_j},\, y_{_k})\in P(\mathring{\Omega}),\\ u_{_{j,k}}^{n}=g^{n}_{_{j,k}} & \qquad \forall (x_{_j},\, y_{_k},n)\in P(\partial\Omega)\times \mathbb{N}_{_0},\\ \end{array}\right. \end{equation}

donde $f_{j,k}=f(x_{_j}, x_{_k})$, $g^{n}_{_{j,k}}=g(x_{_j},y_{_k},t_{_n})$, $P(\mathring{\Omega})=P(\Omega)\cap\mathring{\Omega}$ y $P(\partial \Omega) = P(\Omega)\cap \partial \Omega$. Observe que $P(\mathring{\Omega})$ representa los puntos internos de la malla y $P(\partial \Omega)$ los puntos de la frontera, así que $P(\Omega)=P(\mathring{\Omega})\cup P(\partial \Omega)$.

El objetivo del FTCS es conocer el valor de $u_{_{j,k}}^{n+1}$, por lo tanto al despejarlo de \eqref{eq:ftcs2} se encuentra la siguiente recurrencia:

\begin{equation}\label{eq:moditer1} u_{_{j,k}}^{n+1} \approx \underbrace{\alpha\Delta t\left(\dfrac{\mathrm{id}}{\alpha\Delta t}+\dfrac{s_{_x}-2\mathrm{id}+a_{_x}}{\Delta x^{2}} + \dfrac{s_{_y}-2\mathrm{id}+a_{_y}}{\Delta y^{2}} \ \right)}_{\mathcal{D}}(u_{_{j,k}}^{n}) = \mathcal{D}(u_{_{j,k}}^{n}). \end{equation} El modelo \eqref{eq:moditer1} sugiere el siguiente algoritmo para calcular la temperatura en la placa $\Omega$.

Método implícito BTCS

Si ahora se aplica BTCS, se encuentra en la ecuación de difusión \eqref{eq: difusión} que:

\begin{equation}\label{eq:btcs} \dfrac{u_{_{j,k}}^{n+1} - u_{_{j,k}}^{n}}{\Delta t}\approx \alpha \left( \dfrac{u_{_{j+1,k}}^{n+1}-2u_{_{j,k}}^{n+1}+u_{_{j-1,k}}^{n+1}}{\Delta x^{2}} + \dfrac{u_{_{j,k+1}}^{n+1}-2u_{_{j,k}}^{n+1}+u_{_{j,k-1}}^{n+1}}{\Delta y^{2}} \right), \end{equation} nuevamente usando las funciones sucesor y antecesor pero esta vez en la ecuación \eqref{eq:btcs}, el PVI discretizado toma la forma:

\begin{equation}\label{eq:btcs2} \left\{\begin{array}{cc} \dfrac{(\mathrm{id}-a_{_t})}{\Delta t}(u_{_{j,k}}^{n+1}) \approx \alpha\left(\dfrac{s_{_x}-2\mathrm{id}+a_{_x}}{\Delta x^{2}} + \dfrac{s_{_y}-2\mathrm{id}+a_{_y}}{\Delta x^{2}}\right)(u_{_{j,k}}^{n+1}) & \qquad \forall (x_{_j},\, y_{_k})\in P(\mathring{\Omega}),\\ u_{_{j,k}}^{0}=f_{_{j,k}} & \qquad \forall (x_{_j},\, y_{_k})\in P(\mathring{\Omega}),\\ u_{_{j,k}}^{n}=g_{_{j,k}}^{n} & \qquad \forall (x_{_j},\, y_{_k},\,n)\in P(\partial\Omega)\times \mathbb{N}_{_0},\\ \end{array}\right. \end{equation} Ahora el objetivo del BTCS es conocer el valor de $u_{_{j,k}}^{n}$, por lo tanto se tiene la siguiente recurrencia:

\begin{equation}\label{eq:moditer2} u_{_{j,k}}^{n} = \underbrace{\alpha\Delta t \left(\dfrac{\mathrm{id}}{\alpha\Delta t}-\dfrac{s_{_x}-2\mathrm{id}+a_{_x}}{\Delta x^{2}} - \dfrac{s_{_y}-2\mathrm{id}+a_{_y}}{\Delta y^{2}} \ \right)}_{\mathcal{D}}(u_{_{j,k}}^{n+1}) = \mathcal{D}(u_{_{j,k}}^{n+1}). \end{equation} Para el modelo \eqref{eq:moditer2} se tiene el siguiente algoritmo:

Ejemplo

Se usará el método FTCS para estimar la distribución de temperatura en una placa y su variación en el tiempo. Se considera una placa de dimensiones $8\,cm$ en el eje $x$ y $5\,cm$ en el eje $y$, en este caso, el dominio $\Omega$ es el arreglo rectangular $[0,8]\times [0,\,5]$. Cuando se aplica una partición de orden $l = 30$ y $m = 30$, los pasos para cada eje son $\Delta x=4/15$ y $\Delta y=1/6$. Considere para todos los nodos de $\Omega$ una temperatura inicial de $25\,$°C y las siguientes condiciones de frontera: Sobre la frontera superior e izquierda una temperatura variable linealmente de $100\,$°C a $25\,$°C. Se realiza la estimación de la distribución de calor en la placa para cada punto $(x_{_{j}},\, y_{_{k}})$ para cada tiempo de $0$ a $2000$ s.

Al implementar la ecuación \eqref{eq:moditer1} en Python con las condiciones anteriores se obtienen los siguientes resultados:

Conclusiones

Con el modelo iterativo \eqref{eq:moditer1} el cálculo de la temperatura para cada punto $(x_{_j},\, y_{_k})$ en el tiempo $t_{_{n+1}}$ depende de los valores en el tiempo $t_{_n}$ de los puntos adyacentes en las direcciones de los ejes $x$ y $y$, junto con el valor en $(x_{_j},\, y_{_k})$ como puede apreciarse en la figura 1. En este caso las iteraciones calculan los valores de $(x_{_j}, y_{_k})$ con la información en un tiempo pasado.

En el caso del modelo iterativo \eqref{eq:moditer2} se invierte el proceso con respecto al tiempo, es decir, el cálculo de la temperatura para cada punto $(x_{_j},\, y_{_k})$ en el tiempo $t_{_{n}}$ depende de los valores en el tiempo $t_{_{n+1}}$ de los puntos adyacentes en las direcciones de los ejes $x$ e $y$, junto con el valor en $(x_{_j},\, y_{_k})$ como puede apreciarse en la figura 2. Y en este caso las iteraciones calculan los valores de $(x_{_j}, y_{_k})$ con la información en un tiempo futuro.

Al momento de inicializar los respectivos algoritmos para FCTS y BCTS, se asignan valores de arranque para la temperatura en el punto $(x_{_j},\, y_{_k})$ para todo tiempo $t$. Esto permite usar información del futuro $t_{_{n+1}}$ o del pasado $t_{_{n-1}}$ al iniciar las iteraciones.

Referencias