21.3: Sistema de dos ODEs lineales de primer orden

Es posible abordar directamente numéricamente el sistema de segundo orden de la Sección 21.2 por ejemplo mucho más general y de hecho es la base para la solución numérica de sistemas de ODEs de prácticamente cualquier tipo.

Representación espacial estatal de ODEs escalares de segundo orden

En esta sección, desarrollamos una representación espacial estatal de la ODE canónica de segundo orden. Recordemos que la ODE de interés es de la forma

\(\frac{d^{2} u}{d t^{2}}+2 \zeta \omega_{n} \frac{d u}{d t}+\omega_{n}^{2} u=\frac{1}{m} f(t), \quad 0<t<t_{f}\),

\ begin {alineado}
u (0) &=u_ {0},\\
\ frac {d u} {d t} (0) &=v_ {0}.
\ end {alineado}

Debido a que se trata de una ecuación de segundo orden, necesitamos dos variables para describir completamente el estado del sistema. Escojamos estas variables de estado para ser

\(w_{1}(t)=u(t) \text { and } w_{2}(t)=\frac{d u}{d t}(t)\),

correspondientes al desplazamiento y velocidad, respectivamente. Tenemos la relación trivial entre\(w_1\) y\(w_2\)

\(\frac{d w_{1}}{d t}=\frac{d u}{d t}=w_{2}\).

Además, la ODE de segundo orden gobernante puede ser reescrita en términos\(w_{1}\) y\(w_{2}\) como

\(\frac{d \omega_{2}}{d t}=\frac{d}{d t} \frac{d u}{d t}=\frac{d^{2} u}{d t^{2}}-2 \zeta \omega_{n} \frac{d u}{d t}=-\omega_{n}^{2} u+\frac{1}{m} f=-2 \zeta \omega_{n} w_{2}-\omega_{n}^{2} w_{1}+\frac{1}{m} f\).

Juntos, podemos reescribir la ODE original de segundo orden como un sistema de dos ODEs de primer orden,

\ (\ dfrac {d} {d t}\ left (\ begin {array} {c}
w_ {1}\\
w_ {2}
\ end {array}\ right) =\ left (\ begin {array} {c}
w_ {2}\
-w_ {n} ^ {2} w_ {1} -2\ zeta w_ {n} w_ {2} +\ frac {1} {m} f
\ end {array}\ right)\).

Esta ecuación se puede escribir en forma de matriz

\ [\ frac {d} {d t}\ left (\ begin {array} {l}
w_ {1}\
w_ {2}
\ end {array}\ right) =\ underbrackets {\ left (\ begin {array} {cc}
0 & 1\\
-\ omega_ {n} ^ {2} & -2\ zeta\ omega_ {n}
\ end array}\ derecha)} _ {A}\ izquierda (\ begin {array} {c}
w_ {1}\
w_ {2}
\ end {array}\ right) +\ left (\ begin {array} {c}
0\\
\ frac {1} {m} f
\ end {array}\ derecha)\ tag {21.2}\]

con la condición inicial

\(w_{1}(0)=u_{0} \text { and } w_{2}(0)=w_{n}\)

Si definimos\(w=\left(w_{1} \quad w_{2}\right)^{\mathrm{T}}\) y\(F=\left(\begin{array}{ll}0 & \frac{1}{m} f\end{array}\right)^{\mathrm{T}}\), entonces

\ [\ frac {d w} {d t} =A w+f,\ quad w (t=0) =w_ {0} =\ left (\ begin {array} {l}
u_ {0}\
v_ {0}
\ end {array}\ derecha)\ tag {21.3}\],

resume sucintamente la representación “estado-espacio” de nuestra ODA

Solución por expansión modal

Para resolver esta ecuación, primero encontramos los valores propios de\(A\). Recordamos que los valores propios son las raíces de la ecuación caracteristaica\(p(\lambda ; A)=\) det (\((\lambda I-A)\), donde det se refiere al determinante. (En la práctica real para sistemas grandes los valores propios no se calculan a partir de la ecuación característica.

En nuestro caso 2×2 obtenemos

\ (p (\ lambda; A) =\ nombreoperador {det} (\ lambda I-A) =\ nombreoperador {det}\ left (\ begin {array} {cc}
\ lambda & -1\\
\ omega_ {n} ^ {2} &\ lambda+2\ zeta\ omega_ {n}
\ end {array}\ derecha) =\ lambda^ {2} +2\ zeta\ omega_ {n}\ lambda+\ omega_ {n} ^ {2}\).

Los valores propios, las raíces de la ecuación característica, son así

\(\lambda_{1,2}=-\zeta \omega_{n} \pm \omega_{n} \sqrt{\zeta^{2}-1}\).

Asumiremos en lo sucesivo que el sistema está subamortiguado (es decir,\(\zeta<1\)), en cuyo caso es más conveniente expresar los valores propios como

\(\lambda_{1,2}=-\zeta \omega_{n} \pm i \omega_{n} \sqrt{1-\zeta^{2}} \)

Tenga en cuenta que como el valor propio tiene una parte imaginaria distinta de cero, la solución será oscilatoria y como la parte real es negativa (a la izquierda del plano complejo) la solución es estable. Consideramos ahora los vectores propios.

Hacia ese fin, primero generalizamos nuestro primer discussson de vectores de componentes de valor real al caso de vectores de componentes de valor complejo. A saber, si se nos dan dos vectores\(v \in \mathbb{C}^{m \times 1}\),\(w \in \mathbb{C}^{m \times 1}-v\) y cada uno\(w\) son vectores de columna con entradas\(m\) complejas,\(-\) el producto interno ahora viene dado por

\[\beta=v^{\mathrm{H}_{w}} w=\sum_{j=1}^{m} v_{j}^{*} w_{j} \tag{21.4} \],

donde\(\beta\) es en general complejo,\(H\) significa Hermitiano (transposición compleja) y reemplaza T por transposición, y "denota conjugado complejo\(-80 v_{j}=\operatorname{Real}\left(v_{j}\right)+i \operatorname{Imag}\left(v_{j}\right)\) y\(v_{j}^{*}=\operatorname{Real}\left(v_{j}\right)-\)
\(i \operatorname{Imag}\left(v_{j}\right)\), para\(i=\sqrt{-1}\).

Los diversos conceptos construidos sobre el producto interno cambian de manera similar. Por ejemplo, dos vectores de valor complejo\(v\) y\(w\) son ortogonales si\(v^{\mathrm{H}} w=0\). Lo más importante es que la norma del vector de valor complejo ahora viene dada por

\[\|v\|=\sqrt{v^{\mathbb{H}_{v}}}=\left(\sum_{j=1}^{m} v_{j}^{*} v_{j}\right)^{1 / 2}=\left(\sum_{j=1}^{m}\left|v_{j}\right|^{2}\right)^{1 / 2} \tag{21.5} \],

donde\(|\cdot|\) denota el módulo complejo;\(\left|v_{j}\right|^{2}=v_{j}^{*} v_{j}=\left(\operatorname{Real}\left(v_{j}\right)\right)^{2}+\left(\operatorname{Imag}\left(v_{j}\right)\right)^{2}\). Obsérvese que la definición (21.5) de la norma asegura que\(\|v\|\) es un número real no negativo, como esperaríamos de una longitud.

Para obtener los vectores propios, debemos encontrar una solución a la ecuación

\[(\lambda I-A) X=0 \tag{21.6} \]

para\(\lambda=\lambda_{1}\left(\Rightarrow\right.\) autovector\(\chi^{1} \in \mathbb{C}^{2}\)) y\(\lambda=\lambda_{2}\) (\(\Rightarrow\)vector propio\(\chi^{2} \in \mathbb{C}^{2}\)). Las ecuaciones (21.6) tendrán una solución ya que se\(\lambda\) ha elegido para hacer\((\lambda I-A)\) singular: las columnas de\(\lambda I-A\) son\(x \neq 0\), de las columnas de las\(\lambda I-A\) cuales ylelds el vector cero.

Procediendo con el primer vector propio, escribimos\(\left(\lambda_{1} I-A\right) \chi^{1}=0\) como

\ (\ left (\ begin {array} {cc}
-\ zeta\ omega_ {n} +i\ omega_ {n}\ sqrt {1-\ zeta^ {2}} & -1\\
\ omega_ {n} ^ {2} &\ zeta\ omega_ {n} +i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}
\ end {array}\ derecha)\ left (\ begin {array} {l}
\ chi_ {1} ^ {1}\\
\ chi_ {2} ^ {1}
\ end {array}\ right) =\ left (\ begin {array} {l}
0\\
0
\ end {array}\ right)\)

para obtener (digamos, establecer\(\chi_{1}^{1}=c\)),

\ (\ chi^ {1} =c\ izquierda (\ begin {array} {c}
1\\
\ frac {-\ omega_ {n} ^ {2}} {\ omega_ {n} +i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}}
\ end {array}\ derecha)\).

Ahora elegimos\(c\) lograr\(\left\|x^{1}\right\|=1\), cediendo

\ (x^ {1} =\ frac {1} {\ sqrt {1+\ omega_ {n} ^ {2}}\ left (\ begin {array} {c}
1\\
-\ zeta\ omega_ {n} +i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}
\ end {array}\ derecha)\).

De manera similar obtenemos\(\left(\lambda_{2} I-A\right) \chi^{2}=0\) del segundo vector propio

\ (x^ {2} =\ frac {1} {\ sqrt {1+\ omega_ {n} ^ {2}}\ left (\ begin {array} {c}
1\\
-\ zeta\ omega_ {n} -i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}
\ end {array}\ derecha)\)

que satisface\(\left\|x^{2}\right\|=1\).

Ahora introducimos dos vectores adicionales,\(\psi^{1}\) y\(\psi^{2}\). El vector\(\psi^{1}\) se elige para satisfacer\(\left(\psi^{1}\right)^{\mathrm{H}} \chi^{2}=0\) y\(\left(\psi^{1}\right)^{\mathrm{H}} \chi^{1}=1\), mientras que el vector\(\psi^{2}\) se elige para satisfacer\(\left(\psi^{2}\right)^{\mathrm{H}^{\mathrm{H}}} \chi^{1}=0\) y\(\left(\psi^{2}\right)^{\mathrm{H}} \chi^{2}=\) 1. Encontramos, después de un poco de álgebra,

\ (\ psi^ {1} =\ frac {\ sqrt {1+\ omega_ {n} ^ {2}}} {2 i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}}\ left (\ begin {array} {c}
-\ zeta\ omega_ {n} +i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}\\
-1
\ end {array}\ derecha),\ psi^ {2} =\ frac {\ sqrt {1+\ omega_ {n} ^ {2}}} {-2 i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}}\ left (\ begin {array} {c}
-\ zeta\ omega_ {n} -i\ omega_ {n}\ sqrt {1-\ zeta^ {2}}\\
-1
\ end {array}\ right)\).

Estas elecciones pueden parecer misteriosas, pero en un momento veremos la utilidad de este sistema “bi-ortogonal” de vectores. (Los pasos aquí de hecho corresponden a la “diagonalización” de\(A\).)

Ahora escribimos\(w\) como una combinación lineal de los dos vectores propios, o “modos”

\[\begin{array}{rl}w(t) & =z_{1}(t) \chi^{1}+z_{2}(t) \chi^{2} \\ & =S z(t) \end{array} \tag{21.7}\]

donde

\(S = \left(\chi^{1} \chi^{2}\right)\)

es la\(2 \times 2\) matriz cuya\(j^{\text {th }}\) -columna viene dada por el\(j^{\text {th-eigenvector, } \chi^{j}}\). A continuación insertamos (21.7) en (21.3) para obtener

\[\chi^{1} \frac{d z_{1}}{d t}+\chi^{2} \frac{d z_{2}}{d t}=A\left(\chi^{1} z_{1}+\chi^{2} z_{2}\right)+F \tag{21.8}\]

\[\left(\chi^{1} z_{1}+\chi^{2} z_{2}\right)(t=0)=w_{0} \tag{21.9}\]

Ahora aprovechamos los\(\psi\) vectores.

Primero multiplicamos (21.8) por\(\left(\psi^{1}\right)^{\mathrm{H}}\) y aprovechamos\(\left(\psi^{1}\right)^{\mathrm{H}} \chi^{2}=0,\left(\psi^{1}\right)^{\mathrm{H}} \chi^{1}=1\), y\(A \chi^{j}=\lambda_{j} \chi^{j}\) para obtener

\[\frac{d z_{1}}{d t}=\lambda_{1} z_{1}+\left(\psi^{1}\right)^{\mathrm{H}} F \tag{21.10}\]

si multiplicamos de manera similar (21.9) obtenemos

\[z_1 \left( t=0 \right) = \left( \psi \right)^{\mathrm{H}}w_0 \tag{21.11}\]

El mismo procedimiento pero ahora con\(\left(\psi^{2}\right)^{\mathrm{H}}\) más que\(\left(\psi^{1}\right)^{\mathrm{H}}\) da

\[\frac{d z_{2}}{d t}=\lambda_{2} z_{2}+\left(\psi^{2}\right)^{\mathrm{H}} F \tag{21.12}\];

\[z_{2}(t=0)=\left(\psi^{2}\right)^{\mathrm{H}} w_{0} \tag{21.13}\]

Por lo tanto, observamos que nuestra expansión modal reduce nuestro sistema\(2 \times 2\) ODE acoplado en dos ODE desacopladas.

El hecho de que\(\lambda_{1}\) y\(\lambda_{2}\) sean complejos significa que\(z_{1}\) y también\(z_{2}\) son complejos, lo que podría parecer inconsistente con nuestra ecuación real original (21.3) y solución real\(w(t)\). No obstante, señalamos que\(\lambda_{2}=\lambda_{1}^{*}\) y\(\psi^{2}=\left(\psi^{1}\right)^{*}\) y por lo tanto\(z_{2}=z_{1}^{*}\). Así pues, de (21.7) se desprende que,\(\chi^{2}=\left(\chi^{1}\right)^{*}\) como también,

\(w=z_{1} \chi^{1}+z_{1}^{*}\left(\chi^{1}\right)^{*}\)

y por lo tanto

\(w=2 \operatorname{Real}\left(z_{1} x^{1}\right)\).

Tras la superposición, nuestra solución es realmente real, según lo deseado. el interés aquí en la descomposición modal es como una forma de entender cómo elegir un esquema de ODE para un sistema de dos (posteriores\(n\)) ODEs, y, para el esquema elegido, cómo elegir\(\Delta t\) para la estabilidad.

Aproximación numérica de un sistema de dos ODEs

Biela-Nicolson

La aplicación del esquema Crank-Nicolson a nuestro sistema (21.3) es idéntica a la aplicación del esquema Crank-Nicolson a una ODE escalar. En particular, tomamos directamente el esquema del ejemplo 21.1.8 y reemplazamos por\(\tilde{w}^j \in \mathbb{R}^2\) y\(\tilde{u}^j \in \mathbb{R}\)\(g\) con\(A \tilde{w}^j + F^j\) para obtener

\[\tilde{w}^{j}=\tilde{w}^{j-1}+\frac{\Delta t}{2}\left(A \tilde{w}^{j}+A \tilde{w}^{j-1}\right)+\frac{\Delta t}{2}\left(F^{j}+F^{j-1}\right) \tag{21.14}\].

(Tenga en cuenta si nuestra fuerza\(f\) es constante en el tiempo entonces\(F^{j}=F\).) En general si se desprende de argumentos de consistencia que obtendremos el mismo orden de convergencia que para el problema escalar - si (21.14) es estable. El tema diffeult para los sistemas es la estabilidad: ¿Un esquema particular tendrá un buen esquema de estabilidad será estable? (Esto último es particularmente importante para esquemas explícitos).

Para abordar estas preguntas volvemos a aplicar el análisis modal pero ahora a nuestras ecuaciones discretas (21.14). En particular, escribimos

\[\tilde{w}^{j}=\tilde{z}_{1}^{j} x^{1}+\bar{z}_{2}^{j} x^{2} \tag{21.15}\],

donde\(\chi^{1}\) y\(x^{2}\) son los vectores propios de\(A\) como se derivan en la sección anterior. Ahora insertamos (21.15) en (21.14) y multiplicamos por\(\left(\psi^{1}\right)^{\mathrm{H}}\) y al\(\left(\psi^{2}\right)^{\mathrm{H}}-\) igual que en la sección anterior\(-\) para obtener

\[z_{1}^{j}=z_{1}^{j-1}+\frac{\lambda_{1} \Delta t}{2}\left(\tilde{z}_{1}^{j}+z_{1}^{j-1}\right)+\left(\psi^{1}\right)^{\mathrm{H}} \frac{\Delta t}{2}\left(F^{j}+F^{j-1}\right) \tag{21.16}\],

\[z_{2}^{j}=z_{2}^{j-1}+\frac{\lambda_{2} \Delta t}{2}\left(z_{2}^{j}+z_{2}^{j-1}\right)+\left(\psi^{2}\right)^{\mathrm{H}} \frac{\Delta t}{2}\left(F^{j}+F^{j-1}\right) \tag{21.17}\],

con las condiciones iniciales correspondientes (que no son relevantes para nuestra discusión actual).

Ahora recordamos que para el problema del modelo

\[\dfrac{d u}{d t}=\lambda u+f \tag{21.18}\],

análogo a (21.10), llegamos al esquema Crank-Nicolson

\[\tilde{u}^{j}=\tilde{u}^{j-1}+\frac{\lambda \Delta t}{2}\left(\tilde{u}^{j}+\tilde{w}^{j-1}\right)+\frac{\Delta t}{2}\left(f^{j}+f^{j-1}\right) \tag{21.19}\],

análogo a (21.16). Trabajando hacia atrás, para (21.19) y por lo tanto (21.16) para ser una aproximación estable a (21.18) y por lo tanto (21.10), debemos exigir\(\lambda \Delta t\), y por lo tanto\(\lambda_{1} \Delta t\), residir en el que la ecuación de diferencia (21.16) explotará -y por lo tanto también (21.14) por virtura de (21.15) - si no\(\lambda_{1} \Delta t\) está en el no sombreado región de la Figura 21.12 (a). Por argumentos similares, también\(\lambda_{2} \Delta t\) debe estar en la región no sombreada de la Figura 21.12 (a). En este caso, sabemos que ambos\(\lambda_{1}\) y\(\lambda_{2}\) -para nuestra ecuación particular, es decir, para nuestra matriz particular\(A\) (que determina los valores propios\(\lambda_{1}\),\(\left.\lambda_{2}\right)\) - están en el plano izquierdo, y por lo tanto en la región de estabilidad absoluta de Crank-Nicolson; así Crank-Nicolson es incondicionalmente estable —estable para todos\(\Delta t\) - para nuestra ecuación particular y convergerá\(O\left(\Delta t^{2}\right)\) como as\(\Delta t \rightarrow 0\).

Destacamos que el procedimiento numérico viene dado por (21.14), y no por (21.16), (21.17). El decompositjon modal es sólo para fines de comprensión y análisis - para determinar si un esquema es estable y si es así para qué valores de\(\Delta t\). (Para una\(2 \times 2\) matriz\(A\) la descomposición modal completa es simple. Pero para sistemas más grandes, como consideraremos en la siguiente sección, la descomposición modal completa es muy cara. De ahí que preferimos discretizar directamente la ecuación original, como en (21.14). Esto que\(\Delta t\) en (21.16) y (21.17) son los mismos -ambos se originan en la ecuación (21.14). Lo discutimos más adelante en el contexto de ecuaciones rígidas.

Receta General

Consideramos ahora un sistema general de\(n=2\) ODEs dado por

\[\frac{d w}{d t}=A w+F \tag{21.20}\],

\[w(0)=w_{0} \tag{21.20}\]

donde\(w \in \mathbb{R}^{2}, A \in \mathbb{R}^{2 \times 2}\) (una\(2 \times 2\) matriz)\(F \in \mathbb{R}^{2}\), y\(w 0 \in \mathbb{R}^{2}\). A continuación discretizamos (21.20) por cualquiera de los esquemas desarrollados anteriormente para la ecuación escalar

\(\dfrac{d u}{d t}=g(u, t)\)

simplemente sustituyendo\(w\) por\(u\) y\(Aw + F\) para\(g(u, t)\). Denotaremos el esquema por\(\mathbb{S}\) y la región de estabilidad absoluta asociada por\(\mathcal{R}_{\mathbb{S}}\). Recordemos que (\ mathcal {R} _ {\ mathbb {S}}\) es el subconjunto del plano complejo que contiene todos\(\lambda \Delta t\) para los cuales el esquema\(\mathbb{S}\) aplicado\(g(u, t) = \lambda u\) es absolutamente estable.

Por ejemplo, si aplicamos el esquema Euler Forward\(\mathbb{S}\) obtenemos

\[\tilde{w}^{j}=\tilde{w}^{j-1}+\Delta t\left(A \tilde{w}^{j-1}+F^{j-1}\right) \tag{21.21}\]

mientras que Euler hacia atrás como\(\mathbb{S}\) rendimientos

\[\tilde{w}^{j}=\tilde{w}^{j-1}+\Delta t\left(A \tilde{w}^{j}+F^{j}\right) \tag{21.22}\],

y Crank-Nicolson como\(\mathbb{S}\) da

\[\tilde{w}^{j}=\tilde{w}^{j-1}+\frac{\Delta t}{2}\left(A \tilde{w}^{j}+A \tilde{w}^{j-1}\right)+\frac{\Delta t}{2}\left(F^{j}+F^{j-1}\right) \tag{21.23}\].

Un esquema de múltiples pasos como AB2 como\(\mathbb{S}\) da

\[\tilde{w}^{j}=\tilde{w}^{j-1}+\Delta t\left(\frac{3}{2} A \tilde{w}^{j-1}-\frac{1}{2} A \tilde{w}^{j-2}\right)+\Delta t\left(\frac{3}{2} F^{j-1}-\frac{1}{2} F^{j-2}\right) \tag{21.24}\]

Los diagramas de estabilidad para estos cuatro esquemas\(\mathcal{R}_{\mathbb{S}}\),, están dados por la Figura 21.9, Figura 21.7, Figura 21.12 (a) y Figura 21.11 (b), respectivamente.

A continuación asumimos que podemos calcular los dos valores propios de\(A, \lambda_{1}\), y\(\lambda_{2}\). Un particular\(\Delta t\) conducirá a un esquema estable si y sólo si los dos puntos\(\lambda_{1} \Delta t\) y\(\lambda_{2} \Delta t\) ambos se encuentran dentro\(\mathcal{R}_{\mathbb{S}}\). Si alguno o ambos de los dos puntos\(\lambda_{1} \Delta t\) o se\(\lambda_{2} \Delta t\) encuentran afuera\(\mathcal{R}_{\mathbb{S}}\), entonces debemos disminuir\(\Delta t\) hasta ambos\(\lambda_{1} \Delta t\) y\(\lambda_{2} \Delta t\) mentir dentro\(\mathcal{R}_{\mathbb{S}}\). El paso de tiempo crítico\(\Delta t_{\mathrm{cr}}\),, se define como el mayor\(\Delta t\) para el cual los dos rayos\(\left[0, \lambda_{1} \Delta t\right],\left[0, \lambda_{2} \Delta t\right]\), ambos se encuentran dentro\(\mathcal{R}_{\mathbb{S}} ; \Delta t_{\mathrm{cr}}\) dependerá de la forma y tamaño de\(\mathcal{R}_{\mathbb{S}}\) y la “orientación” de los dos rayos\(\left[0, \lambda_{1} \Delta t\right],\left[0, \lambda_{2} \Delta t\right]\).

Podemos derivar\(\Delta t_{\mathrm{cr}}\) de una manera ligeramente diferente. Primero definimos\(\widehat{\Delta t}_{1}\) que es el más grande de\(\Delta t\) tal manera que\(\left[0, \lambda_{1} \Delta t\right]\) está el rayo\(\mathcal{R}_{\mathbb{S}}\);\(\widehat{\Delta t}_{2}\) a continuación definimos que sea el más grande de\(\Delta t\) tal manera que el rayo\(\left[0, \lambda_{2} \Delta t\right]\) esté adentro\(\mathcal{R}_{\mathbb{S}}\). Entonces podemos deducir eso\(\Delta t_{\mathrm{cr}}=\min \left(\widehat{\Delta t}_{1}, \widehat{\Delta t}_{2}\right)\). En particular, observamos que si\(\Delta t>\Delta t_{\mathrm{cr}}\) entonces uno de los dos modos -y de ahí toda la solución- explotará. También podemos ver aquí de nuevo la dificultad con ecuaciones rígidas en las que\(\lambda_{1}\) y\(\lambda_{2}\) son muy diferentes:\(\widehat{\Delta t} t_{1}\) pueden ser (digamos) mucho más grandes que\(\widehat{\Delta t}_{2}\), pero\(\widehat{\Delta t}_{2}\) dictarán\(\Delta t\) y así nos obligarán a dar muchos pasos de tiempo -muchos más de los requeridos para resolver los más lentos (menor\(\left|\lambda_{1}\right|\) asociado con decaimiento más lento u oscilación más lenta) que suele ser el comportamiento de interés.

En lo anterior asumimos, como es casi siempre el caso, que los\(\lambda\) están en el plano izquierdo. Para cualquier $\ lambda$ que se encuentre en el plano derecho, nuestra condición es volteada: ahora debemos asegurarnos de que los no\(\lambda \Delta t\) están en la región de estabilidad absoluta para obtener las soluciones de crecimiento (inestables) deseadas.

Cerremos esta sección con dos ejemplos.

Ejemplo 21.3.1 Sistema de masa de resorte sin amortiguar

En este ejemplo, revisamos el sistema de masa de resorte sin amortiguar considerado en la sección anterior. Los dos valores propios de\(A\) son\(\lambda_{1}=i \omega_{n}\) y\(\lambda_{2}=i \omega_{n}\); sin pérdida de generalidad, nos fijamos\(\omega_{n}=1.0\). Consideraremos la aplicación de varios esquemas de integración numérica diferentes al problema; para cada integrador, evaluamos su aplicabilidad con base en la teoría (apelando al diagrama de estabilidad absoluta) y verificaremos nuestra evaluación a través de experimentos numéricos.

\((i)\)Euler Forward es una mala elección ya que ambos\(\lambda_{1} \Delta t\) y\(\lambda_{2} \Delta t\) están afuera\(\mathcal{R}_{\mathbb{S}=E F}\) para todos\(\Delta t\). El resultado del experimento numérico, mostrado en la Figura 21.20 (a), confirma que la amplitud de la oscilación crece para ambos\(\Delta t=0.5\) y\(\Delta t=0.025\); el paso de tiempo menor da como resultado una amplificación más pequeña (artificial).

\((ii)\)Euler Hacia atrás también es una mala elección ya que $\ lambda_ {1}\ Delta t$ y $\ lambda_ {2}\ Delta t$ están en el interior de $\ mathcal {R} _ {\ mathrm {S} =\ text {EB}} $ para todos $\ Delta t$ y por lo tanto la solución discreta decaerá a pesar de que la solución exacta sea una oscilación no en descomposición. La figura 21.20 (b) confirma la valoración.

\((iii)\)Crank-Nicolson es una muy buena opción ya que\(\lambda_{1} \Delta t \in \mathcal{R}_{\mathrm{S}=\mathrm{CN}}, \lambda_{2} \Delta t \in \mathcal{R}_{\mathrm{S}=\mathrm{CN}}\) para todos\(\Delta t\), y además se\(\lambda_{1} \Delta t, \lambda_{2} \Delta t\) encuentran en el límite de\(\mathcal{R}_{S=\mathrm{CN}}\) y por lo tanto la solución discreta, al igual que la solución exacta, no se desintegrará. La Figura 21.20 (c) confirma que Crank-Nicolson conserva los resultados en un error de fase notable.

\((iv)\)Cuatro etapas Runge-Kutta (RK4) es una opción razonablemente buena desde entonces\(\lambda_{1} \Delta t\) y se\(\lambda_{2} \Delta t\) encuentran cerca del límite de\(\mathcal{R}_{\mathbb{S} = \rm RK4}\) for\(| \lambda_i \Delta t | \lesssim 1\). La Figura 21.20 (d) muestra que, para el problema considerado, RK4 sobresale no solo en preservar la amplitud de la oscilación sino también en alcanzar la fase correcta.

Obsérvese en el análisis anterior que el diagrama de estabilidad absoluta sirve no solo para determinar la estabilidad sino también la naturaleza de la solución discreta en cuanto al crecimiento, o decaimiento, o incluso estabilidad neutra, sin crecimiento ni decaimiento. (Esto último no implica que la solución discreta sea exacta, ya que además de demostrar en particular la presencia de errores de fase en ausencia de errores de amplitud).

Figura 21.20: Comparación de esquemas de integración numérica para un sistema de masa de resorte sin amortiguar con\(ω_n = 1.0\).

Ejemplo 21.3.2 Sistema amortiguador-masa-amortiguador sobreamortiguado: un sistema rígido de ODEs

En nuestro segundo ejemplo, consideramos un sistema de resorte-masa-amortiguador (muy) sobreamortiguado con\(\omega_{n}=1.0\) y\(\zeta=100\). Los valores propios asociados con el sistema son

\ begin {alineado}
&\ lambda_ {1} =-\ zeta\ omega_ {n} +\ omega_ {n}\ sqrt {\ zeta^ {2} -1} =-0.01\\
&\ lambda_ {2} =-\ zeta\ omega_ {n} -\ omega_ {n}\ sqrt {\ zeta^ {2} -1} =-99.99.
\ end {alineado}

Como antes, perturbaremos el sistema por un desplazamiento inicial unitario. El modo lento con $\ lambda_ {1} =-0.01$ dicta la respuesta del sistema. Sin embargo, para esquemas condicionalmente estables, la estabilidad se rige por el modo rápido con\(\lambda_2 = −99.99\). De nuevo consideramos cuatro integradores de tiempo diferentes: dos explícitos y dos implícitos.

\((i)\)Euler Forward es estable para\(\Delta t \lesssim 0.02\) (es decir\(\Delta t_{\mathrm{cc}}=2 /\left|\lambda_{2}\right|\)). La Figura 21.21 (a) muestra que el esquema rastrea con precisión la solución exacta (bastante benigna) para\(\Delta t=0.02\), pero se vuelve inestable y diverge exponencialmente para\(\Delta t=0.0201\). Así, el paso de tiempo máximo es limitado (dictado por\(\lambda_{2}\)). En otras palabras, aunque la respuesta del sistema es benigna, no podemos usar grandes pasos tíme para ahorrar en costos computacionales.

\((ii)\)Similar al caso Euler Forward, el esquema de cuatro etapas Runge-Kutta (RK4) exhibe un comportamiento exponencialmente divergente para\(\Delta t>\Delta t_{\mathrm{cr}} \approx 0.028\), como se muestra en la Figura 21.21 (b). El paso de tiempo máximo está nuevamente limitado por la estabilidad.

\((iii)\)Euler Hacia atrás es incondicionalmente estable, y así la elección del paso de tiempo viene dictada por la capacidad de aproximar la respuesta del sistema, que viene dictada por\(\lambda_{1}\). Figura 21.21 (c) grande\(\Delta t=5.0\) ya que la respuesta del sistema es bastante lenta.

\((iv)\)Crank-Nicolson también es incondicionalmente estable. Para el mismo conjunto de pasos de tiempo, Crank-Nicolson produce una aproximación más precisa que Euler Hacia atrás, como se muestra en la Figura 21.21 (d), debido a su precisión de orden superior.

En la comparación anterior, los esquemas incondicionalmente estables requirieron muchos menos pasos de tiempo (y bence mucho menos esfuerzo computacional) que los esquemas condicionalmente estables. Por ejemplo, Crank Nicolson con\(\Delta t=5.0\) requiere aproximadamente 200 veces menos pasos de tiempo que el esquema RK4 (con una elección estable del paso de tiempo). Más importante aún, como la escala de tiempo más corta (es decir, el valor propio más grande) dicta la estabilidad, los esquemas condicionalmente estables no permiten al usuario usar pasos de tiempo grandes incluso si los modos rápidos no son de interés para el usuario. Como se mencionó anteriormente, los sistemas rígidos son omnipresentes en la ingeniería, y los ingenieros a menudo no están interesados en la escala de tiempo más pequeña presente en el sistema. (Recordemos el ejemplo de la escala de tiempo asociada a la dinámica de un jet de pasajeros y la asociada a remolinos turbulentos; a menudo los ingenieros solo están interesados en caracterizar la dinámica de la aeronave, no los remolinos). En estas situaciones, los esquemas incondicionalmente estables permiten a los usuarios elegir un paso de tiempo apropiado independientemente de las limitaciones de estabilidad.

___________________________ . ___________________________

Para terminar, queda claro incluso a partir de estos sencillos ejemplos que un esquema explícito de propósito general idealmente incluiría alguna parte tanto del eje real negativo como del eje imaginario. Los esquemas que exhiben este comportamiento incluyen AB3 y RK4. De estos dos esquemas, a menudo se prefiere RK4 debido a una gran región de estabilidad; también RK4, un método multietapa, no sufre los problemas de puesta en marcha que a veces complican las técnicas de varios pasos.