§15 Método de MaCormack & Paullay (1972)

A continuación explicaremos brevemente un método numérico para resolver temporalmente la ecuación:

$$\frac{ \partial r }{ \partial t } = A \frac{ \partial r }{ \partial \theta } + B,$$ (15.40)

para un ángulo polar $\theta$ fijo, donde en general $A$ y $B$ son funciones que dependen de la posición $r$, del ángulo $\theta$ y del tiempo $t$. Adjunta con la ec.(3.40) se impone la condición inicial:

$$r \bigg\vert _{ t = 0 } = r_0,$$ (15.41)

para todo $\theta$. Generalizando el método de Euler (Abramowitz & Stegun, 1970), la solución avanzada de la ec.(3.40) es:

$$r \left( \theta, t + \Delta t \right) = r ( \theta, t ) + \frac{ \partial r }{ \partial t } \bigg\vert _{ \theta, t }\Delta t,$$ (15.42)

donde $\Delta t$ es un pequeño intervalo de tiempo. Con ayuda de la ec.(3.40) la ec.(3.42) puede escribirse como:

$$r \left( \theta, t + \Delta t \right) = r ( \theta, t ) + \left( A \frac{ \partial r }{ \partial \theta } + B \right)_{ \theta, t } \Delta t.$$ (15.43)

En la práctica sucede que el cálculo de $\partial r / {\partial \theta}$ debe hacerse de manera discreta, que no se considera en la ec.(3.43). Un método que considera el cálculo de la derivada de la posición con respecto del ángulo polar discretamente es el método de MaCormack & Paullay (1972), que describiremos a continuación.

El problema en cuestión está caracterizado por las cantidades $A, B, r, \theta$. De aquí se pueden construir dos tiempos: $r / \left\vert B \right\vert$ y $\Delta \theta / \left\vert A\right\vert$, donde $\Delta \theta = \left( \theta_\text{max} - \theta_\text{min} \right) / n$ y $\theta_$max$$ y $\theta_$min$$ se refieren a los valores del ángulo polar máximo y mínimo respectivamente, $n$ es un número entero suficientemente grande. De esta manera, el intervalo de tiempo $\Delta t$ se elige como:

$$\Delta t = \zeta \left\{ \frac{ r ( \theta_i ) }{ \left\vert B \right\vert }, \frac{ \Delta \theta }{ \left\vert A \right\vert } \right\}_{ i = 0 }^{ i = n},$$ (15.44)

donde $\zeta$ es un número suficientemente pequeño y $\theta_i = i \Delta \theta + \theta_$min$$ ( $i = 0, 1, \ldots n$).

Para calcular el valor de la posición $r$ al avanzar el tiempo $\Delta t$, primero se utiliza la ec.(3.43) calculando la derivada de manera discreta hacia la derecha, excepto en el borde, donde se calcula hacia la izquierda, obteniendo así el valor $r_$d$( \theta_i )$. A esta posición avanzada se le avanza nuevamente en el intervalo de tiempo $\Delta t$, utilizando la ec.(3.43) y calculando la derivada de manera discreta hacia la izquierda, excepto en el borde, donde se calcula hacia la derecha, obteniendo el valor $r_$l$( \theta_i )$. Se puede mostrar que la posición:

$$r( \theta_i, t + \Delta t ) = \frac{ 1 }{ 2 } \left\{ r_\text{d} ( \theta_i ) + r_\text{l} ( \theta_i ) \right\},$$ (15.45)

coincide con la solución de la ec.(3.40) expandida en serie de potencias hasta segundo orden temporal y espacialmente.

En la práctica suele suceder que la solución obtenida con la ec.(3.45) presenta picos. Para suavizar estas soluciones se utiliza un método conocido como corrector de flujo (Boris & Book, 1973; Boris & Book, 1976), que describiremos a continuación. Consideremos:

$$R( \theta_i, t + \Delta t ) = r( \theta_i, t + \Delta t ) + \frac... ...\theta_{ i + 1 }, t + \Delta t ) - r( \theta_{ i - 1 }, t + \Delta t )\right\},$$

$$R( \theta_0, t + \Delta t ) = \frac{ 1 }{ 2 } \left\{ r( \theta_0, t + \Delta t ) + r( \theta_1, t + \Delta t ) \right\},$$

$$R( \theta_n, t + \Delta t ) = \frac{ 1}{ 2 } \left\{ r( \theta_{ n - 1 }, t + \Delta t ) + r( \theta_{n}, t + \Delta t ) \right\}.$$ (15.46)

Para analizar la existencia de picos en la solución de la ec.(3.45), hagamos $q_i \equiv r( \theta_i, t + \Delta t ) - r( \theta_{ i - 1 }, t + \Delta t )$ y analicemos distintos casos:

1. $q_i q_{ i + 1 } < 0$ (existe pico). Se escoge la solución suavizada y avanzada en el tiempo de la ec.(3.40) como:

   $$S( \theta_i, t + \Delta t ) = R( \theta_i, t + \Delta t ),$$ (15.47)

   
   
   con $i=0,1,2,...,n$.
2. $q_i q_{i+1} \geq 0$ (no existe pico). Se escoge la solución poco suavizada o no suavizada de la ec.(3.40) como:

   $$S( \theta_i, t + \Delta t ) = R( \theta_i, t + \Delta t ) - \frac... ...theta_{ i + 1 }, t + \Delta t ) - r( \theta_{ i - 1 }, t + \Delta t ) \right\},$$ (15.48)

   
   
   donde $\epsilon = 1$ es una solución no suavizada y $\epsilon \lesssim 1$ es una solución poco suavizada para $i \not = 0, n$ (nosotros utilizamos el valor $\epsilon = 0.9$). En los bordes ($i = 0, n$) la solución se toma como en la ec.(3.46).