13.5: Simulación de Modelos de Campo Continuos
- Page ID
- 115865
La simulación de modelos de campo continuos escritos en PDEs no es una tarea fácil, ya que implica fácilmente una enorme cantidad de cómputos si se desea obtener resultados de simulación bastante precisos, y también porque ciertos tipos de dinámicas espaciales a veces pueden causar comportamientos altamente sensibles. De hecho, la mayoría de las “supercomputadoras” construidas para computación científica 2 se utilizan para resolver problemas de PDE a gran escala, por ejemplo, para simular sistemas climáticos globales, propiedades complejas de materiales, etc.
Cubrir todas las técnicas avanzadas de simulaciones numéricas de PDEs está mucho más allá del alcance de este libro de texto introductorio. En cambio, intentaremos un enfoque mucho más sencillo. Como ya discutimos los autómatas celulares (CA), convertiremos modelos basados en PDE en CA discretizando el espacio y el tiempo, para luego simularlos simplemente como modelos CA con variables de estado continuas.
Discretizar el tiempo en una PDE no es diferente de lo que discutimos en la Sección 6.3. Al reemplazar la derivada del tiempo por su análogo discretizado, la PDE original
\[\dfrac{\partial{f}}{\partial{t}}= F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}}} , \dots, x, t)\label{(13.29)} \]
se convierte
\[\dfrac{\Delta{f}}{\Delta{t}} =\dfrac{f(x, t +\Delta{t}) -f(x, t)}{\Delta{t}}=F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}}} , \dots, x, t), \label{(13.30)} \]
\[f(x, t+ \Delta{t}) =f(x,t) +F(f, \dfrac{\partial{f}}{\partial{x}}, \dfrac{\partial^{2}{f}}{\partial{x^{2}} }, \dots, x, t)\Delta{t} \label{(13.31)} \]
que ahora es una ecuación de diferencia a lo largo del tiempo (¡fíjese que arriba no es un laplaciano!).
Aún no hemos terminado, porque el lado derecho aún puede contener derivados espaciales (\(∂f/∂x\),\(∂^2f/∂x^2\),...). También necesitamos discretizarlos sobre el espacio, para poder desarrollar una versión CA de esta PDE. Aquí hay una forma de discretizar una derivada espacial de primer orden:
\[\dfrac{\partial{f}}{\partial{x}} \approx \dfrac{\Delta{f}}{\Delta{x}} =\dfrac{f(x, t +\Delta{t}) -f(x, t)}{\Delta{x}} \label{(13.32)} \]
Esto no es incorrecto matemáticamente, pero hay una cuestión práctica. Esta discretización hace que el espacio sea asimétrico a lo largo del eje x, porque la derivada espacial\(x\) depende de\(f(x + ∆x,t)\), pero no de\(f(x−∆x,t)\). A diferencia del tiempo que tiene una clara dirección asimétrica del pasado al futuro, el eje espacial debería modelarse mejor simétricamente entre la izquierda y la derecha (o arriba y abajo). Por lo tanto, una mejor manera de realizar la discretización espacial es la siguiente:
\[\dfrac{\partial{f}}{\partial{x}} \approx \dfrac{2\Delta{f}}{2\Delta{x}}= \dfrac{f(x, t +\Delta{t}) -f(x, t)}{2\Delta{x}}\label{(13.33)} \]
Esta versión trata simétricamente a la izquierda y a la derecha.
Del mismo modo, podemos derivar la versión discretizada de una derivada espacial de segundo orden, de la siguiente manera:
\[ \dfrac{\partial^{2}{f}}{\partial{x^{2}}} \approx \dfrac{\dfrac{\partial{f}}{\partial{x}}|_{x+\Delta{x}} -\dfrac{\partial{f}}{\partial{x}}|_{x-\Delta{x}}}{2\Delta{x}} \nonumber \]
\[f(x + ∆x + ∆x,t ) \qquad f(x−∆x + ∆x,t) \nonumber \]
\[ \approx \dfrac{\dfrac{-f(x +\Delta{x} -\Delta{x,t})}{ 2\Delta{x}} - \dfrac{-f(x-\Delta{x} -\Delta{x,t})}{2\Delta{x}}} {2\Delta{x}} \nonumber \]
\[ =\dfrac{f(x+2\Delta{x,t}) +f(x-2\Delta{x,t}) -2f(x,t)}{2\Delta{x^{2}}} \nonumber \]
\[=\dfrac{f(x+\Delta{x, t}) +f(x-\Delta{x, t}) -2f(x,t)}{\Delta{x^{2}}} \nonumber \]
\[ =\dfrac{f(x +\Delta{x, t}) +f(x-\Delta{x, t}) -2f(x,t)}{\Delta{x^{2}}} \text{by replacing 2∆x → ∆x } \label{(13.34)} \]
Además, podemos usar este resultado para discretizar a un laplaciano de la siguiente manera:
\[\nabla^{2} f =\dfrac{\partial^{2} {f}}{\partial{x}^{2}_{1}} + \dfrac{\partial^{2} {f}}{\partial{x}^{2}_{2}} +.....+ \dfrac{\partial^{2} {f}}{\partial{x}^{2}_{n}} \nonumber \]
\[\approx \dfrac{f(x_1 + ∆x_1,x_2,...,x_n,t) + f(x_1 −∆x_1,x_2,...,x_n,t)−2f(x_1,x_2,...,x_n,t) }{∆x^2_1 } \nonumber \]
\[ + \dfrac{f(∆x_1,x_2,...,x_n,t) + f(∆x_1,x_2 -\Delta{x_2},...,x_n,t)−2f(x_1,x_2,...,x_n,t) }{∆x^2 _2} \nonumber \]
\[+\dots \nonumber \]
\[ + f(x_1,x_2,...x_n +\Delta{x_n},t) + f(x_1,x_2,...,x_n -\Delta{x_n},t) \nonumber \]
\[= ( f(x_1+ \Delta{x}, x_{2},\dots, x_{n}, t ) +f(x_1 -\Delta{x}, x_2, ...x_n, t) -2f(x_{1}, x_{2},...., x_{n}, t) \nonumber \]
\[+ f(x_1, x_2, +\Delta{x}, ...., x_n , t)+f(x_1, x_2, -\Delta{x}, ...,x_n, t) \nonumber \]
\[+.... \nonumber \]
\[ + f(x_1, x_2, ...., x_n +\Delta{x, t}) +f(x_1, x_2,..., x_n -\Delta{x, t}) \nonumber \]
\[-2nf(x_1, x_2,...,x_n, t)/ {\Delta{x^{2}} } \label{(13.35)} \]
\[ \text {(by replacing ∆xi → ∆x for all i)} \nonumber \]
El numerador de la fórmula anterior tiene una interpretación intuitiva: “Simplemente agregue los valores de\(f\) en todos sus vecinos más cercanos (por ejemplo, arriba, abajo, izquierda y derecha para un espacio 2-D) y luego restar su propio valor\(f(x,t)\) tantas veces como el número de esos vecinos”. O, de manera equivalente: “Mide la diferencia entre tu prójimo y tú mismo, y luego resuma todas esas diferencias”. Es interesante ver que un operador de orden superior como los Lapacianos puede aproximarse mediante un conjunto tan simple de operaciones aritméticas una vez que el espacio está discretizado.
También se pueden derivar derivadas espaciales de orden superior de manera similar, pero no creo que necesitemos cubrir esos casos aquí. En fin, ahora tenemos un conjunto básico de herramientas matemáticas (Eqs. \ ref {(13.31)},\ ref {(13.33)},\ ref {(13.34)},\ ref {(13.35)}) para discretizar modelos basados en PDE para que podamos simular su comportamiento como modelos de CA.
Trabajemos en un ejemplo. Considera discretizar una ecuación de transporte simple sin fuente o sumidero en un espacio 2-D, donde la velocidad de las partículas viene dada por un vector constante:
\[ \dfrac{\partial{c}}{\partial{t}} =-\nabla \cdot(cw) \label{(13.36)} \]
\[ w = \binom{w_x}{w_y} \label{(13.37)} \]
Podemos discretizar esta ecuación tanto en el tiempo como en el espacio, de la siguiente manera:
\[c(x, y, t +\delta{t}) \approx c(x, y,t) -\nabla \cdot (c(x,y,t)w)\delta{t} \label{(13.38)} \]
\[=c(x, y, t ) - {\begin {pmatrix} \dfrac{\partial}{\partial{x}} \\ \dfrac{\partial}{\partial{y}} \end{pmatrix}}^{T} (c(x,y,t) \begin{pmatrix} w_x \\ w_y \end{pmatrix})\Delta{t} \label{(13.39)} \]
\[=c(x, y, t ) - (w_x \dfrac{\partial{c}}{\partial{x}} + w_y \dfrac{\partial{c}}{\partial{y}}) \Delta{ t} \label{(13.40)} \]
\[\approx c(x, y, t ) - (w_x \dfrac{c(x +\Delta{h, y, t}) - c(x-\Delta{h,y, t})}{2\Delta{h}} + w_y \dfrac{c(x,y +\Delta{h, t}) - c(x, y-\Delta{h, t})}{2\Delta{h}}) \Delta{t} \label{(13.41)} \]
\[= c(x, y, t ) -(w_{x}c(x +\Delta{h, y, t}) - w_{x}c(x -\Delta{h, y, t}) + w_{y}c(x,y +\Delta{h, t}) - w_{y}c(x, y-\Delta{h, t})) \dfrac{\Delta{t}}{2\Delta{h}} \label{(13.42)} \]
Esto ahora está completamente discretizado tanto para el tiempo como para el espacio, y por lo tanto se puede usar como una función de transición de estado de CA. Podemos implementar este modelo de CA en Python, usando Code 11.5 como plantilla. Aquí hay un ejemplo de implementación:
En este ejemplo, simulamos la ecuación de transporte con\((w_x,w_y) = (0.01,0.03)\) en un espacio 2-D [0,1] × [0,1] con condiciones de contorno periódicas, partiendo de una configuración inicial que tiene un pico en el medio del espacio:
\[c(x,y, 0) =exp(-\dfrac{(x-0.5)^{2} +(y+ 0.5)^{2}}{0.2^{2}}) \label{(13.43)} \]
Las resoluciones temporales y espaciales se establecen en\(∆t = 0.01\) y\(∆h = 0.01\), respectivamente (de manera que el tamaño de la cuadrícula sea 100×100). Este es un proceso de simulación muy lento, por lo que el StepSize se establece en 50 para omitir visualizaciones innecesarias (ver la última línea). La gráfica de superficie se utiliza para visualizar la configuración en ejes tridimensionales. Tenga en cuenta que es necesario llamar al comando show al final de la función observe, ya que los cambios realizados en los ejes 3D no se reflejan automáticamente a la visualización (al menos en la implementación actual de matplotlib).
El resultado de esta simulación se muestra en la Fig. 13.5.1 donde se puede ver claramente el pico siendo transportado a una dirección dada por\((w_x,w_y)\). Además, gracias a la naturaleza interactiva de pycxsimulator.py, esta visualización 3-D ahora es manipulable de forma interactiva incluso si no la estás ejecutando desde Anaconda Spyder, Enthought Canopy u otros entornos interactivos. Puede hacer clic y arrastrar sobre él para rotar la superficie y observar su estructura 3-D desde varios ángulos.
Si no necesitas una parcela de superficie 3-D, yo ucan usa imshow en su lugar. Esto hace que el código sea un poco más simple y, sin embargo, el resultado puede ser tan bueno como antes (Fig. 13.5.2):
Una cosa de la que debo advertirles es que el método de simulación discutido anteriormente todavía se basa en el método forward de Euler, por lo que acumula fácilmente errores numéricos. Este es el mismo tema de la estabilidad de las soluciones y la posibilidad de “artefactos” que discutimos en la Sección 6.4, que ahora surgen de la discretización tanto del espacio como del tiempo. En el ejemplo de la ecuación de transporte anterior, tales problemas pueden ocurrir si elige un valor que es demasiado grande para\(∆t\) o\(∆h\), relativo a la escala espacial/temporal del comportamiento del sistema (que está determinada por w en este caso). Por ejemplo, si aumentas\(∆t\) a 0.1, obtendrás el
resultado mostrado en la Fig. 13.5.3. Esto puede parecer genial, pero en realidad es un resultado no válido. Se puede decir que no es válido si entiende la naturaleza de la ecuación de transporte; solo debe transportar la superficie y no cambiar su forma. Pero aunque no estés seguro de si el resultado que obtuviste es válido o no, puedes intentar incrementar las resoluciones temporales/espaciales (es decir, hacer\(∆t\) y\(∆h\) más pequeñas) para ver si el resultado varía. Si el resultado no se ve afectado por esto, entonces usted sabe que las resoluciones espaciales/temporales originales ya eran suficientes y que el resultado original era así válido. Pero si no, el resultado original probablemente contenía errores numéricos, por lo que es necesario seguir aumentando las resoluciones hasta que el resultado se vuelva consistente.
Antes de pasar al siguiente tema, me gustaría señalar que los métodos de discretización/simulación discutidos anteriormente pueden extenderse a PDEs que involucran múltiples variables de estado. Este es un ejemplo de este tipo: interacciones entre dos poblaciones distribuidas espacialmente, un escape de la difusión\(b\):
\[\dfrac{\partial{a}}{\partial{t}} =- \mu_{a} \nabla \cdot(a(-\nabla{b})) \label{(13.44)} \]
\[ \dfrac{\partial{b}}{\partial{t}} =\mu_{b} \nabla^{2} b \label{(13.45)} \]
Aquí\(µ_a\) y\(µ_b\) están los parámetros que determinan la movilidad de las dos especies. Estas ecuaciones se pueden expandir y luego discretizar para un espacio 2-D, de la siguiente manera (¡también deberías intentar hacerlo todo por ti mismo!) :
\[\dfrac{\partial{a}}{\partial{t}} =-\mu_{a} (\begin{pmatrix} \dfrac{\partial}{\partial{x}} \\ \dfrac{\partial}{\partial{y}}) \end{pmatrix}^{T} \begin{pmatrix} -a\dfrac{\partial{b}}{\partial{x}} \\ -a\dfrac{∂b}{∂y}\end{pmatrix} \label{(13.46)} \]
\[=\mu_{a} (\dfrac{\partial}{\partial{x}} (a\dfrac{∂b}{∂x}) + \dfrac{\partial}{\partial{y}}(a\dfrac{∂b}{∂y})) \label{(13.47)} \]
\[= \mu_{a} ( \dfrac{\partial{a}}{\partial{x}} \dfrac{∂b}{∂x} + a\dfrac{∂^{2}b}{∂x^{2}} + \dfrac{\partial{a}}{\partial{y}} \dfrac{∂b}{∂y} +a\dfrac{∂^{2}b}{∂y^{2}} \label{(13.48)} \]
\[=\mu_{a} (\dfrac{\partial{a}}{\partial{x}} \dfrac{∂b}{∂x} +\dfrac{\partial{a}}{\partial{y}} \dfrac{∂b}{∂y} a\nabla^{2}b) \\label{(13.49)} \]
\ [a (x, y, t +\ Delta {t})\ aprox a (x, y, t) + µ_a (\ dfrac {a (x + ∆h, y, t) −a (x−∆h, y, t)} {2∆h}
\ dfrac {b (x + ∆h, y, t) −b (x−∆h, y, t)} {2∆h} +\ dfrac {a (x, y + ∆h, t) −a (x, y−∆h, t)} {2∆h}
\ dfrac {b (x, y + ∆h, t) −b (x, y−∆h, t)} {2∆h}
+ a (x, y, t) ×\ dfrac {b (x + ∆h h, y, t) + b (x−∆h, y, t) + b (x, y + ∆h, t) + b (x, y−∆h, t) −4b (x, y, t)} {\ Delta {h^ {2}}})\ Delta {t}\ etiqueta {(13.50)}\]
\[a'c \approx a_c +\mu_{a} (\dfrac{a_R −a_L}{2∆h} \dfrac{b_R −b_L}{2∆h} + \dfrac{a_V −a_D}{2∆h} \dfrac{b_V −b_D}{2∆h} + a_C \dfrac{b_R + b_L + b_U + b_D −4b_C} {∆h^{2}} )∆t \label{(13.51)} \]
Obsérvese que utilicé notaciones más simples en la última ecuación, donde los subíndices\((C, R, L, U, D)\) representan estados de la celda central así como sus cuatro vecinos, y\(a'C\) es el siguiente valor de\(a_C\). Mientras tanto, la ecuación discretizada para b viene dada simplemente por
\[b'_C \approx b_C +\mu_b \dfrac{b_R + b_L + b_U + b_D −4b_C}{\Delta{h^{2}}} \Delta{t}. \label{(13.52)} \]
A continuación se muestra un código de ejemplo para simular este modelo PDE comenzando con una configuración inicial con dos picos, uno para a y otro para\(b\):
Tenga en cuenta que la función observe ahora usa subplot para mostrar dos campos escalares para\(a\) y\(b\). El resultado se muestra en la Fig. 13.5.4.
Modificar el Código 13.7 para que ambas especies se difundan y traten de escapar unas de otras. Ejecuta simulaciones para ver cómo se comporta el sistema.
2 Si no tienes claro de lo que estoy hablando aquí, consulta http://top500.org/.