Análisis de estabilidad de Von Neumann

La estabilidad de los esquemas numéricos está estrechamente relacionada con el error numérico. Un esquema de diferencias finitas es estable si los errores cometidos en un paso temporal del cálculo no provocan que los errores se amplifiquen a medida que continúan los cálculos. Un «esquema neutro estable» es aquel en el que los errores permanecen constantes a medida que avanzan los cálculos. Si los errores disminuyen y finalmente se amortiguan, se dice que el esquema numérico es estable. Si, por el contrario, los errores aumentan con el tiempo, se dice que el esquema numérico es inestable. La estabilidad de los esquemas numéricos puede investigarse 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 sea acotada. En general, la estabilidad puede ser difícil de investigar, especialmente cuando la ecuación considerada es no lineal.

En ciertos casos, la estabilidad de von Neumann es necesaria y suficiente para la estabilidad en el sentido de Lax-Richtmyer, tal y como se utiliza en el teorema de equivalencia de Lax. La ecuación diferencial parcial y los modelos de esquema de diferencias finitas son lineales; la ecuación diferencial parcial es de coeficiente constante con condiciones de frontera periódicas y solo tiene dos variables independientes; y el esquema no utiliza 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 utiliza en lugar de un análisis de estabilidad más detallado para proporcionar una buena estimación de las restricciones (si las hay) sobre los tamaños de paso utilizados en el esquema debido a su relativa simplicidad.

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

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

$${\displaystyle u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)}$$

(1)

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

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

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

$${\displaystyle \epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)}$$

(2)

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 decaimiento 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 ${\displaystyle x}$, en el intervalo ${\displaystyle L}$, como

$${\displaystyle \epsilon (x,t)=\sum _{m=-M}^{M}E_{m}(t)e^{{i}k_{m}x}}$$

(3)

donde el número de onda ${\displaystyle k_{m}={\frac {\pi m}{L}}}$ con ${\displaystyle m=-M,\dots ,-2,-1,0,1,2,\dots ,M}$ y ${\displaystyle M=L/\Delta x}$. La dependencia temporal del error se incluye suponiendo que la amplitud del error ${\displaystyle E_{m}}$ 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.

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

$${\displaystyle \epsilon (x,t)=\int _{-{\frac {\pi }{\Delta x}}}^{\frac {\pi }{\Delta x}}E_{k}(t)e^{ikx}dk}$$

(4)

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

$${\displaystyle \epsilon _{m}(x,t)=E_{m}(t)e^{ik_{m}x}}$$

(5a)

si se utiliza una serie de Fourier o

$${\displaystyle \epsilon _{k}(x,t)=E_{k}(t)e^{ikx}}$$

(5b)

si se utiliza una integral de Fourier.

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

Las características de estabilidad pueden estudiarse 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 observar que $${\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}}}$$ para obtener (tras la simplificación)

$${\displaystyle {\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}=1+r\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right).}$$

(6)

Introducimos ${\displaystyle \theta =k_{m}\Delta x\in [-\pi ,\pi ]}$ y utilizamos las identidades $${\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}}}$$ La ecuación (6) se puede escribir como

$${\displaystyle {\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}=1-4r\sin ^{2}(\theta /2)}$$

(7)

Defina el factor de amplificación

$${\displaystyle G\equiv {\frac {E_{m}(t+\Delta t)}{E_{m}(t)}}}$$

(8)

La condición necesaria y suficiente para que el error permanezca acotado es que ${\displaystyle |G|\leq 1.}$

Por lo tanto, a partir de las ecuaciones (7) y (8), la condición para la estabilidad viene dada por

$${\displaystyle \left|1-4r\sin ^{2}(\theta /2)\right|\leq 1}$$

(9)

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

$${\displaystyle 4r\sin ^{2}(\theta /2)\leq 2}$$

(10)

Para que la condición anterior se cumpla para todos los ${\displaystyle m}$ (y, por lo tanto, para todos los ${\displaystyle \sin ^{2}(\theta /2)}$). El valor más alto que puede tomar el término sinusoidal es 1 y, para esa elección concreta, 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 r={\frac {\alpha \Delta t}{\left(\Delta x\right)^{2}}}\leq {\frac {1}{2}}}$$

(11)

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

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