Estabilidad numérica de Euler Forward para una ecuación diferencial.

Question

Estoy tratando de controlar que la siguiente differentialequation es estable para el método numérico de integración temporal: Euler Forward con $\Delta t = 1$ .

$$ y''(t)+y'(t)+\sin(y(t)) = 0, \ \ \ \ y(0 ) = 1, y'(1) = 0. $$

Utilizando lo siguiente, podemos transformar esto en un sistema de ecuaciones diferenciales. Como sabemos que nuestra EDO es no lineal, utilizaremos la matriz de Jacobi de $\bf{y'}= \begin{bmatrix} y_1'(t) \\ y_2'(t) \end{bmatrix}$ . \begin{align} y_1&=y &\Rightarrow y_1' &= y_2,\\ y_2&=y' &\Rightarrow y_2'& = -y_2-\sin(y_1). \end{align} La matriz de Jacobi es:

\begin{equation} J(y_1,y_2)= \begin{bmatrix} 0 & 1\\ -\cos(y_1) & -1 \end{bmatrix} . \fin{s} de la ecuación

Vamos a calcular los valores propios de esta matriz, estos son: \begin{align} \lambda_1 &= \frac{-1 + \sqrt{1-4 \cos(y_1)} }{2},\\ \lambda_2 &= \frac{-1 - \sqrt{1-4 \cos(y_1)} }{2}. \end{align}

Ahora, mi pregunta es: Estos valores propios pueden ser complejos, entonces cómo es posible que Euler Forward sea estable, recuerdo que Euler Forward no es estable si los valores propios tienen partes complejas.

Gracias por su tiempo,

K. Kamal

Answer 1

La condición de estabilidad es que para los valores propios $λ$ con parte real negativa tienes $$|1+hλ|<1.$$ Para su problema estos valores propios son $λ=\frac12(-1\pm i\sqrt{3})$ y $λ=-\frac12(1+\sqrt5)$ para que la condición sea $$ (1-\tfrac12h)^2+\tfrac34h^2<1\iff h^2-h<0\text{ and }h(1+\sqrt5)<4 $$ y como $h>0$ el límite de estabilidad es $h<\min(1, \frac{4}{1+\sqrt5})=1$ .

Sin embargo, esto sólo significa que se evita un comportamiento visiblemente erróneo de la solución numérica, es decir, que si la solución exacta se mueve hacia un punto fijo, la solución numérica se aleje de él. El error global de tamaño $O(h)$ o por tiempo variable, $O((e^{Lt}-1)h)$ sigue en pie. Puede ser aconsejable tomar $h$ mucho menor que el límite de estabilidad para obtener un error global razonable.

Mientras que la solución exacta converge muy rápidamente hacia el origen, con un poco más de una rotación visible alrededor de cero, la solución numérica con tamaño de paso $h=0.5$ todavía muestra un poco esa rapidez mientras $h=0.75$ tiene una convergencia bastante lenta y $h=1$ parece indeciso si se acerca al origen en absoluto.