Análisis de estabilidad de von Neumann

Procedimiento de análisis numérico

En análisis numérico , el análisis de estabilidad de von Neumann (también conocido como análisis de estabilidad de Fourier) es un procedimiento utilizado para verificar la estabilidad de los esquemas de diferencias finitas aplicados a ecuaciones diferenciales parciales lineales . ^[1] El análisis se basa en la descomposición de Fourier del error numérico y fue desarrollado en el Laboratorio Nacional de Los Álamos después de haber sido descrito brevemente en un artículo de 1947 por los investigadores británicos John Crank y Phyllis Nicolson . ^[2] Este método es un ejemplo de integración temporal explícita donde la función que define la ecuación gobernante se evalúa en el momento actual. Más tarde, el método recibió un tratamiento más riguroso en un artículo ^[3] coescrito por John von Neumann .

Estabilidad numérica

La estabilidad de los esquemas numéricos está estrechamente asociada con el error numérico . Un esquema de diferencias finitas es estable si los errores cometidos en un paso de tiempo del cálculo no hacen que los errores se magnifiquen a medida que continúan los cálculos. Un esquema neutralmente estable es uno en el que los errores permanecen constantes a medida que se realizan los cálculos. Si los errores decaen y finalmente desaparecen, se dice que el esquema numérico es estable. Si, por el contrario, los errores crecen con el tiempo, se dice que el esquema numérico es inestable. La estabilidad de los esquemas numéricos se puede investigar realizando un análisis de estabilidad de von Neumann. Para problemas dependientes del tiempo, la estabilidad garantiza que el método numérico produzca una solución acotada siempre que la solución de la ecuación diferencial exacta esté acotada. La estabilidad, en general, puede ser difícil de investigar, especialmente cuando la ecuación en consideración es no lineal .

En ciertos casos, la estabilidad de von Neumann es necesaria y suficiente para la estabilidad en el sentido de Lax-Richtmyer (como se usa en el teorema de equivalencia de Lax ): los modelos de EDP y de esquema de diferencias finitas son lineales; la EDP es de coeficientes constantes con condiciones de contorno periódicas y tiene solo dos variables independientes; y el esquema no usa más de dos niveles de tiempo. ^[4] La estabilidad de von Neumann es necesaria en una variedad mucho más amplia de casos. A menudo se usa en lugar de un análisis de estabilidad más detallado para proporcionar una buena estimación de las restricciones (si las hay) en los tamaños de paso utilizados en el esquema debido a su relativa simplicidad.

Ilustración del método

El método de von Neumann se basa en la descomposición de los errores en series de Fourier . Para ilustrar el procedimiento, considere la ecuación de calor unidimensional definida en el intervalo espacial , con la notación donde son los valores x específicos y son la secuencia de valores t . $${\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}$$ ${\displaystyle L}$ $${\displaystyle u_{j}^{n}=u(x_{j},t^{n})}$$ ${\displaystyle x_{j}}$ ${\displaystyle t^{n}}$

Podemos discretizar la ecuación de calor ^[5] como

dónde $${\displaystyle r={\frac {\alpha \,\Delta t}{\left(\Delta x\right)^{2}}}}$$

Entonces, la solución de la ecuación discreta se aproxima a la solución analítica de la EDP en la cuadrícula. ${\displaystyle u_{j}^{n}}$ ${\displaystyle u(x,t)}$

Defina el error de redondeo como donde es la solución de la ecuación discretizada ( 1 ) que se calcularía en ausencia de error de redondeo, y es la solución numérica obtenida en aritmética de precisión finita . Dado que la solución exacta debe satisfacer exactamente la ecuación discretizada, el error también debe satisfacer la ecuación discretizada. ^[6] Aquí asumimos que también satisface la ecuación (esto solo es cierto en precisión de máquina). Por lo tanto ${\displaystyle \epsilon _{j}^{n}}$ $${\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}$$ ${\displaystyle u_{j}^{n}}$ ${\displaystyle N_{j}^{n}}$ ${\displaystyle u_{j}^{n}}$ ${\displaystyle \epsilon _{j}^{n}}$ ${\displaystyle N_{j}^{n}}$

es una relación de recurrencia para el error. Las ecuaciones ( 1 ) y ( 2 ) muestran que tanto el error como la solución numérica tienen el mismo comportamiento de crecimiento o decrecimiento con respecto al tiempo. Para ecuaciones diferenciales lineales con condición de contorno periódica, la variación espacial del error puede expandirse en una serie de Fourier finita con respecto a , en el intervalo , como ${\displaystyle x}$ ${\displaystyle L}$

donde el número de onda con y . La dependencia temporal del error se incluye suponiendo que la amplitud del error es una función del tiempo. A menudo se supone que el error crece o decae exponencialmente con el tiempo, pero esto no es necesario para el análisis de estabilidad. ${\displaystyle k_{m}={\frac {\pi m}{L}}}$ ${\displaystyle m=-M,\dots ,-2,-1,0,1,2,\dots ,M}$ ${\displaystyle M=L/\Delta x}$ ${\displaystyle E_{m}}$

Si la condición de contorno no es periódica, entonces podemos utilizar la integral de Fourier finita con respecto a : ${\displaystyle x}$

Como la ecuación diferencial del error es lineal (el comportamiento de cada término de la serie es el mismo que el de la serie misma), basta considerar el crecimiento del error de un término típico:

Si se utiliza una serie de Fourier o

si se utiliza una integral de Fourier.

Como la serie de Fourier puede considerarse un caso especial de la integral de Fourier, continuaremos el desarrollo utilizando las expresiones para la integral de Fourier.

Las características de estabilidad se pueden estudiar utilizando solo esta forma para el error sin pérdida de generalidad. Para averiguar cómo varía el error en pasos de tiempo, sustituya la ecuación ( 5b ) en la ecuación ( 2 ), después de notar que para obtener (después de la simplificación) $${\displaystyle {\begin{aligned}\epsilon _{j}^{n}&=E_{m}(t)e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=E_{m}(t+\Delta t)e^{ik_{m}x}\\\epsilon _{j+1}^{n}&=E_{m}(t)e^{ik_{m}(x+\Delta x)}\\\epsilon _{j-1}^{n}&=E_{m}(t)e^{ik_{m}(x-\Delta x)},\end{aligned}}}$$

Introduciendo y utilizando las identidades, la ecuación ( 6 ) puede escribirse como ${\displaystyle \theta =k_{m}\Delta x\in [-\pi ,\pi ]}$ $${\displaystyle \sin \left({\frac {\theta }{2}}\right)={\frac {e^{i\theta /2}-e^{-i\theta /2}}{2i}}\qquad \rightarrow \qquad \sin ^{2}\left({\frac {\theta }{2}}\right)=-{\frac {e^{i\theta }+e^{-i\theta }-2}{4}}}$$

Definir el factor de amplificación

La condición necesaria y suficiente para que el error permanezca acotado es que Así, de las ecuaciones ( 7 ) y ( 8 ), la condición de estabilidad viene dada por ${\displaystyle |G|\leq 1.}$

Nótese que el término siempre es positivo. Por lo tanto, para satisfacer la ecuación ( 9 ): ${\displaystyle 4r\sin ^{2}(\theta /2)}$

Para que la condición anterior se cumpla para todos (y, por lo tanto, para todos los ), el valor más alto que puede tomar el término sinusoidal es 1 y, para esa elección en particular, si se cumple la condición del umbral superior, también se cumplirá para todos los puntos de la cuadrícula, por lo que tenemos ${\displaystyle m}$ ${\displaystyle \sin ^{2}(\theta /2)}$

La ecuación ( 11 ) proporciona el requisito de estabilidad para el esquema FTCS aplicado a la ecuación de calor unidimensional. Indica que para un determinado , el valor permitido de debe ser lo suficientemente pequeño para satisfacer la ecuación ( 10 ). ${\displaystyle \Delta x}$ ${\displaystyle \Delta t}$

Un análisis similar muestra que un esquema FTCS para advección lineal es incondicionalmente inestable.

Referencias

1. Análisis de métodos numéricos por E. Isaacson, HB Keller
2. Crank, J.; Nicolson, P. (1947), "Un método práctico para la evaluación numérica de soluciones de ecuaciones diferenciales parciales de tipo de conducción de calor", Proc. Camb. Phil. Soc. , 43 : 50–67, doi :10.1007/BF02127704
3. Charney, JG; Fjørtoft, R.; von Neumann, J. (1950), "Integración numérica de la ecuación de vorticidad barotrópica", Tellus , 2 (4): 237–254, doi : 10.3402/tellusa.v2i4.8607
4. Smith, GD (1985), Solución numérica de ecuaciones diferenciales parciales: métodos de diferencias finitas, 3.ª ed. , págs. 67–68
5. en este caso, utilizando el esquema de discretización FTCS
6. Anderson, JD Jr. (1994). Dinámica de fluidos computacional: conceptos básicos con aplicaciones . McGraw Hill .