Ecuación rígida

En matemáticas computacionales , una ecuación rígida es un problema de valor inicial.

${\displaystyle {\dot {u}}=f(u)\,,\qquad u(0)=u_{0}\,,\qquad t\in [0,T]\,,}$

dónde${\displaystyle f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$, que requiere métodos de paso de tiempo implícitos específicos para su integración numérica eficiente. La caracterización matemática más simple de una ecuación rígida es la condición necesaria

${\displaystyle {\frac {T}{d}}{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\ll -1\,.}$

Desde${\displaystyle {\big (}{\mathrm {div} }_{u}\,f{\big )}(u)={\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)={\mathrm {trace} }\,f'(u)}$, dónde${\displaystyle f'(u)\in {\mathbb {R} }^{d\times d}}$es la matriz jacobiana de${\displaystyle f}$en ese punto${\displaystyle u}$El criterio anterior se evalúa fácilmente y cuantifica la rigidez. A continuación, se deriva, explica e ilustra dicho criterio para ecuaciones rígidas no lineales.

Para un sistema lineal con coeficientes constantes ${\displaystyle {\dot {u}}=Au}$La divergencia es constante, lo que convierte la rigidez en una característica global cuya magnitud está relacionada con la escala de tiempo.${\displaystyle T}$.

En un sistema no lineal , la rigidez suele variar en el espacio y el tiempo a lo largo de la trayectoria de la solución.${\displaystyle u(t)}$donde el criterio cuantifica la rigidez localmente. En cálculos prácticos, las ecuaciones rígidas se resuelven invariablemente utilizando métodos adaptativos. ^{[ 1 ] [ 2 ]}

Existe una amplia bibliografía sobre ecuaciones diferenciales rígidas, pero las descripciones intuitivas y las heurísticas son mucho más comunes que los intentos de una definición rigurosa del concepto. Hairer y Wanner ^{[ 3 ]} describen concisamente la característica más obvia:

Las ecuaciones rígidas son problemas para los que los métodos explícitos no funcionan.

Esto se refiere a la observación de que los métodos de integración explícita se ven obligados a utilizar pasos de tiempo extremadamente pequeños.${\displaystyle h}$para mantener la estabilidad numérica, evitando que dichos métodos sean competitivos. Aunque cada paso es económico, el número total de pasos${\displaystyle N=T/h}$se vuelve prohibitivamente grande y la integración más allá de${\displaystyle [0,T]}$efectivamente se detiene.

En cambio, los métodos implícitos para ecuaciones rígidas requieren una costosa resolución algebraica de ecuaciones en cada paso. El trabajo adicional por paso se compensa con propiedades de estabilidad superiores, lo que permite el uso de pasos de tiempo grandes. Sin restricciones de estabilidad limitantes, el esfuerzo computacional total es manejable y se puede lograr la precisión requerida sin pérdida de eficiencia. Para algunas ecuaciones rígidas, la eficiencia puede ser varios órdenes de magnitud mayor que la de incluso los mejores métodos explícitos.

El fenómeno de rigidez es análogo a problemas de rendimiento bien conocidos en otras áreas del análisis numérico. Por ejemplo, en optimización , los métodos de gradiente tienen limitaciones similares. Utilizar el método de descenso más pronunciado para minimizar una función convexa${\displaystyle F(x)}$conduce a la iteración

${\displaystyle x^{k+1}=x^{k}-h\,{\mathrm {grad} }_{x}F(x^{k})\,,}$

dónde${\displaystyle h}$es el tamaño del paso de búsqueda lineal. Este es simplemente el método de Euler explícito para la ecuación diferencial.${\displaystyle {\dot {x}}=-{\mathrm {grad} }_{x}F(x^{k})}$. Los problemas con gradientes pronunciados son conocidos por su lenta convergencia, debido a las restricciones en la elección de${\displaystyle h}$La solución habitual consiste en sustituir el método del gradiente por algún método de tipo Newton, que esencialmente corresponde al uso de un método "implícito" para ecuaciones diferenciales.

Otro ejemplo se encuentra en la resolución de ecuaciones no lineales de la forma${\displaystyle x=g(x)}$Iteración simple de punto fijo ,${\displaystyle x^{k+1}=g(x^{k})}$, no converge para problemas donde la constante de Lipschitz${\displaystyle L[g]}$es grande, pero solo para contracciones. Sin embargo, la ecuación puede tener una solución única bajo condiciones mucho más débiles pero comunes, por ejemplo, cuando el mapa${\displaystyle gI}$es monótono y la constante de Lipschitz logarítmica de${\displaystyle g}$Satisface${\displaystyle M[g]<1}$Dado que las iteraciones de punto fijo divergen, es necesario utilizar una iteración de tipo Newton.

Un tercer ejemplo son los solucionadores iterativos para la ecuación de Poisson (discretizada).${\displaystyle -\Delta u=f}$Aquí, el método de Jacobi (que intenta evitar el álgebra matricial) es equivalente a utilizar el método explícito de Euler (en pseudotiempo) para la discretización correspondiente del método de líneas de la ecuación de difusión.${\displaystyle u_{t}=\Delta u+f}$, buscando su solución estacionaria. La ecuación de difusión es un ejemplo prototípico de una ecuación rígida: el paso de tiempo de la recursión explícita está severamente restringido por una condición CFL y debe ser proporcional al cuadrado del ancho de la malla en el espacio. Esto hace que la convergencia de la iteración de Jacobi sea extremadamente lenta. Como antes, la solución es usar alguna forma de resolución de ecuaciones algebraicas, por ejemplo, usando un precondicionador o un método multigrid .

Así pues, existen numerosos problemas matemáticos en los que la eficiencia computacional exige el uso de métodos numéricos específicos, generalmente de considerable complejidad. Las ecuaciones rígidas constituyen un ejemplo particular de los problemas de valor inicial en ecuaciones diferenciales ordinarias. No obstante, gracias a un solucionador adaptativo para ecuaciones rígidas bien diseñado, la integración numérica eficiente de la mayoría de estas ecuaciones se ha convertido en una tarea rutinaria.

La primera mención de ecuaciones rígidas se encuentra en Curtiss y Hirschfelder en 1952 ^{[ 4 ]} . Los autores discuten una ecuación de modelo escalar equivalente a

${\displaystyle {\dot {x}}={\frac {1}{a(t,x)}}{\big (}xg(t){\big )}}$

y caracterizar la rigidez de la siguiente manera:

Si${\displaystyle \Delta t}$es la resolución deseada de${\displaystyle t}$o el intervalo que se utilizará en la integración numérica, la ecuación es "rígida" si

${\displaystyle \left|{\frac {a(t,x)}{\Delta t}}\right|\ll 1}$

y${\displaystyle g}$se comporta bien.

Esta caracterización relaciona una escala de tiempo.${\displaystyle \Delta t}$a la tasa de decaimiento de los transitorios, regida por el coeficiente${\displaystyle a(t,x)}$, que se supone negativo. Como${\displaystyle a(t,x)}$Varía en el tiempo y el espacio, lo que permite que la rigidez varíe en consecuencia. Prothero y Robinson introdujeron un problema modelo similar pero más instructivo, ^{[ 5 ]} tomando esencialmente${\displaystyle a(t,x)=1/\lambda <0}$constante y estudiando la ecuación

${\displaystyle {\dot {x}}=\lambda {\big (}xg(t){\big )}+{\dot {g}}(t)\,,\qquad x(0)=x_{0}\neq g(0)\,,\qquad t\in [0,T].}$

La solución consiste en un transitorio decreciente${\displaystyle {\big (}x_{0}-g(0){\big )}\,{\mathrm {e} }^{t\lambda }}$junto con la solución particular${\displaystyle g(t)}$, que se supone que es acotada y suave. La solución matemática${\displaystyle x(t)}$Entonces, obviamente, está acotado. Al contrastar los métodos de Euler explícitos e implícitos, se puede investigar bajo qué circunstancias los métodos numéricos replican las propiedades de la solución matemática. Se considerará el criterio de Curtiss-Hirschfelder para${\displaystyle \Delta t=T}$y${\displaystyle \Delta t=h}$, correspondientes a las condiciones${\displaystyle T\lambda \ll -1}$y${\displaystyle h\lambda \ll -1}$, respectivamente.

El método explícito de Euler para una ecuación diferencial${\displaystyle {\dot {u}}=f(u)}$se define por

${\displaystyle u_{n+1}=u_{n}+hf(u_{n})\,,}$

mientras que el método implícito de Euler se define por la recursión.

${\displaystyle u_{n+1}=u_{n}+hf(u_{n+1})\,.}$

Aquí${\displaystyle h}$es el paso de tiempo (constante), y${\displaystyle u_{n}}$se aproxima a la solución exacta${\displaystyle u(t_{n})}$en ese momento${\displaystyle t_{n}=nh}$. En el método explícito${\displaystyle u_{n+1}}$se calcula directamente mediante una evaluación del campo vectorial.${\displaystyle f}$En el método implícito, sin embargo, hay que resolver la ecuación "algebraica".${\displaystyle u_{n+1}-hf(u_{n+1})=u_{n}}$para${\displaystyle u_{n+1}}$Esto mejora la estabilidad de la recursión.

Para el problema de Prothero-Robinson, el método explícito de Euler produce la solución numérica (en forma integrada).

${\displaystyle x_{n}=\left(1+h\lambda \right)^{n}x_{0}+\sum _{k=0}^{n-1}\left(1+h\lambda \right)^{n-k-1}h{\big (}{\dot {g}}(t_{k})-\lambda g(t_{k}){\big )}\,.}$

El impedimento es claramente visible: potencias positivas del factor${\displaystyle 1+h\lambda }$producir crecimiento exponencial (inestabilidad numérica) a menos que impongamos la condición de estabilidad.${\displaystyle |1+h\lambda |\leq 1}$. Desde${\displaystyle \lambda }$es real y negativo, esto implica${\displaystyle -2\leq h\lambda \leq 0}$Por lo tanto, el método explícito de Euler no puede utilizar tamaños de paso adaptados a la regularidad de${\displaystyle g(t)}$. En cambio, el número total de pasos para completar la integración,${\displaystyle N=T/h\sim T|\lambda |}$, se vuelve extremadamente grande. La eficiencia es inversamente proporcional a la rigidez cuantificada por el producto${\displaystyle T|\lambda |\gg 1}$.

En consecuencia, el método explícito de Euler puede tardar "una eternidad" en resolver un problema escalar con una solución matemática perfectamente suave.${\displaystyle \sim g(t)}$una vez que el transitorio (rápido) se ha atenuado. Esto también ocurre incluso si el valor inicial es${\displaystyle x(0)=g(0)}$Por lo tanto, la restricción del tamaño del paso no es una cuestión de si el transitorio está presente . Se debe al hecho de que cualquier tamaño de paso${\displaystyle h}$de tal manera que${\displaystyle |1+h\lambda |>1}$inevitablemente provoca inestabilidad numérica .

Si en cambio se utiliza el método implícito de Euler, la solución numérica se convierte en

${\displaystyle x_{n}=\left(1-h\lambda \right)^{-n}x_{0}+\sum _{k=0}^{n-1}\left(1-h\lambda \right)^{k-n}h{\big (}{\dot {g}}(t_{k+1})-\lambda g(t_{k+1}){\big )}\,.}$

Aquí, el requisito de estabilidad es${\displaystyle |1-h\lambda |^{-1}\leq 1}$, lo cual se satisface para todos${\displaystyle h>0}$cuando${\displaystyle {\mathrm {Re} }\,\lambda <0}$, sin importar cuán grande${\displaystyle h|\lambda |}$es. Un problema menor es que en la práctica uno necesita resolver el transitorio (esto puede requerir inicialmente un tamaño de paso pequeño), pero una vez que se ha atenuado, se pueden usar tamaños de paso "grandes", en el sentido de que${\displaystyle h}$se adapta a la regularidad de${\displaystyle g(t)}$, independientemente de la magnitud de${\displaystyle |h\lambda |}$De los dos métodos anteriores, solo el método implícito de Euler puede manejar ecuaciones rígidas de manera eficiente.

Ejemplo Consideremos el problema de la prueba de Prothero-Robinson

${\displaystyle {\dot {x}}\,=\,\lambda (x-\sin \omega t)+\omega \cos \omega t\,;\qquad x(0)=x_{0}\,,}$

con solución exacta${\displaystyle x(t)=x_{0}\,{\mathrm {e} }^{\lambda t}+\sin \omega t}$. Tomando${\displaystyle x_{0}=0}$La solución homogénea no está presente en la solución exacta.${\displaystyle x(t)=\sin \omega t}$, pero el término exponencial aparecerá en las soluciones locales que pasan por los puntos${\displaystyle x_{n}}$generado por los métodos numéricos. Tomando${\displaystyle \omega =\pi }$El problema se resuelve en${\displaystyle [0,1]}$, usando${\displaystyle N=5}$y${\displaystyle 10}$pasos, respectivamente (${\displaystyle h=0.2}$y${\displaystyle h=0.1}$). El experimento demuestra el impacto de la variación${\displaystyle h\lambda }$, lo que afecta la estabilidad/amortiguación de los métodos numéricos. Tomamos${\displaystyle \lambda =-20}$, utilizando la misma configuración computacional para los métodos de Euler explícitos e implícitos.

Los resultados se muestran en la Figura 1. Con${\displaystyle T\lambda =-20}$, hay una amortiguación exponencial moderadamente fuerte, como se evidencia en las soluciones locales de rápida contracción. La configuración se ha simplificado deliberadamente con fines de visualización; en ecuaciones rígidas reales, los valores típicos serían mucho mayores. Sin embargo,${\displaystyle \lambda }$es suficientemente grande y negativo como para demostrar las limitaciones del tamaño de paso del método explícito de Euler. Para${\displaystyle N=5}$, su solución calculada se vuelve oscilatoria y diverge de la solución exacta, aunque el valor inicial se tomó en la solución exacta. Repitiendo el mismo ejercicio para${\displaystyle T\lambda =-2000}$Por ejemplo, muestra cómo la diferencia entre los dos métodos aumenta con la rigidez.

En problemas más realistas, la inestabilidad suele ser más dramática y se la denomina "explosión" de la solución. Aquí, el tamaño máximo de paso estable (${\displaystyle h=0.1}$) solo ha sido superado por un pequeño factor. Mientras tanto, el método implícito de Euler procede sin pérdida de estabilidad y es capaz de producir resultados precisos incluso para el tamaño de paso mayor.

En cálculos reales, las ecuaciones rígidas se resuelven mediante métodos adaptativos. Estos seleccionan automáticamente el tamaño del paso para controlar el proceso computacional, generalmente regulando tanto la precisión como la estabilidad. Por lo tanto, si se utiliza un método explícito adaptativo, la solución puede no volverse inestable, ya que la adaptabilidad toma medidas para no exceder el tamaño máximo de paso estable. Sin embargo, este método tiene el inconveniente (prohibitivo) de una baja eficiencia.

La esencia del análisis es que un parámetro del problema,${\displaystyle \lambda }$, que es característico del campo vectorial y está relacionado con el comportamiento transitorio de las soluciones, está vinculado a una escala de tiempo, ya sea el tamaño del paso.${\displaystyle h}$o el rango de integración${\displaystyle T}$directamente. Los productos${\displaystyle h\lambda }$o${\displaystyle T\lambda }$Son adimensionales, invariantes de escala y negativas. Un concepto general de rigidez debe reflejar estas propiedades y, además, ser capaz de abordar la rigidez en ecuaciones diferenciales no lineales.

El problema escalar de Prothero-Robinson puede extenderse a sistemas de ecuaciones lineales,

${\displaystyle {\dot {x}}=A{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad A\in {\mathbb {R} }^{d\times d}}$

así como sistemas no lineales,

${\displaystyle {\dot {x}}=f{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$

siempre que${\displaystyle f(0)=0}$El concepto general de rigidez necesita especificar qué se entiende por un campo vectorial "grande y negativo". El análisis anterior también muestra que la función suave${\displaystyle g}$como mucho desempeña un papel secundario. Por lo tanto, al sustituir${\displaystyle u}$para${\displaystyle x-g(t)}$Tenemos tres tipos de problemas,

${\displaystyle {\begin{aligned}{\dot {u}}\,&=\,\lambda u\\{\dot {u}}\,&=\,Au\\{\dot {u}}\,&=\,f(u)\,.\end{aligned}}}$

Aquí es necesario caracterizar la rigidez de las matrices.${\displaystyle A}$así como para mapas no lineales${\displaystyle f}$La primera ecuación anterior es la conocida ecuación de prueba lineal de Dahlquist , utilizada principalmente para determinar la región de estabilidad de un método de avance temporal. La región de estabilidad${\displaystyle S}$de un método es el conjunto de todos${\displaystyle h\lambda \in \mathbb {C} }$De tal manera que el método numérico produce soluciones acotadas. La discusión previa sobre la rigidez se basa, en efecto, en la ecuación de la prueba de Dahlquist.

Para considerar el segundo sistema de ecuaciones lineales, es común utilizar los valores propios.${\displaystyle \lambda _{k}}$de la matriz${\displaystyle A}$Después de todo, para que un método numérico produzca soluciones acotadas (estables)${\displaystyle {\dot {u}}=Au}$, es necesario y suficiente que todos los autovalores satisfagan${\displaystyle h\lambda _{k}\in S}$. La discusión anterior continúa considerando un valor propio a la vez. El problema lineal es entonces rígido si todos los valores propios tienen partes reales negativas y al menos uno de ellos satisface${\displaystyle T\lambda _{k}\ll -1}$; este último valor propio provocará entonces las restricciones típicas del tamaño del paso que experimenta un método explícito, debido a que todos los métodos explícitos tienen regiones de estabilidad limitadas y requieren${\displaystyle h\lambda _{k}\in S}$. El tamaño del paso debe seleccionarse lo suficientemente pequeño para cumplir con este requisito, véase la Figura 2 .

Sin embargo, es común encontrar una caracterización diferente de la rigidez, en términos de la "relación de rigidez", definida por

${\displaystyle {\frac {\max _{k}|{\mathrm {Re} }\,\lambda _{k}|}{\min _{k}|{\mathrm {Re} }\,\lambda _{k}|}}\,,}$

suponiendo que${\displaystyle \lambda _{k}\in {\mathbb {C} }^{-}}$. ^{[ 6 ]} Aunque los sistemas rígidos suelen tener una gran relación de rigidez, Byrne y Hindmarsh enfatizan que esto no es ni necesario ni suficiente para${\displaystyle {\dot {u}}=Au}$ser una ecuación rígida. ^{[ 7 ]}

En primer lugar, no hay comparación con la escala de tiempo.${\displaystyle T}$, pero solo una comparación intrínseca de "constantes de tiempo"${\displaystyle 1/\lambda _{k}}$donde una gran relación de rigidez simplemente implica que hay componentes transitorios "más rápidos" y "más lentos", a menudo denominados escalas de tiempo muy diferentes. En segundo lugar, como se mencionó anteriormente, existen ecuaciones rígidas escalares; estas tienen una relación de rigidez${\displaystyle 1}$. Y en tercer lugar, la relación de rigidez deja de ser válida para problemas con una matriz singular.${\displaystyle A}$, lo cual es común en algunas aplicaciones. (Un ejemplo es la cinética de reacciones químicas, donde la singularidad corresponde a la conservación de la masa).

Por lo tanto, la relación de rigidez es engañosa: no logra caracterizar los aspectos más importantes de las ecuaciones rígidas .

En la misma línea, se ha sugerido que el problema no lineal es rígido si la matriz jacobiana${\displaystyle f'(u)}$tiene una gran relación de rigidez. ^{[ 8 ]} Naturalmente, esto también falla, como señalan Artemiev y Averina: ^{[ 9 ]}

Por ejemplo, un sistema de EDO autónomo no lineal con${\displaystyle f(0)=0}$se puede decir que es rígido si las partes reales [sic] de los valores propios de la matriz de Jacobi [${\displaystyle f'(0)}$] satisfacen las condiciones anteriores. Sin embargo, la famosa ecuación de van der Pol, que se usa frecuentemente como ejemplo de EDO rígidas, no satisface esta definición... es imposible determinar la rigidez solo mediante los valores propios de la matriz jacobiana para sistemas no lineales.

Aunque parcialmente válido, este argumento también es erróneo, ya que la solución cero de la ecuación de van der Pol es un equilibrio inestable . Por lo tanto, la rigidez no se produce cerca del origen, sino solo a lo largo del ciclo límite, como se verá más adelante, donde se analiza la ecuación de van der Pol con todo detalle. En consecuencia, una caracterización precisa de la rigidez debe ser local, capturando el comportamiento de perturbaciones (pequeñas) cerca de la solución real.

Sustituyendo la relación de rigidez, el enfoque moderno consiste en construir un elemento funcional.${\displaystyle s_{2}[\cdot ]}$, llamado indicador de rigidez , de tal manera que los criterios${\displaystyle s_{2}[A]\!\cdot \!T\ll -1}$y${\displaystyle s_{2}[f]\!\cdot \!T\ll -1}$, respectivamente, cuantifican la rigidez relacionando una propiedad distinta del campo vectorial con la escala de tiempo.${\displaystyle T}$Esta caracterización se basa en observaciones adicionales. Así, Dekker y Verwer ^{[ 10 ]} escriben (énfasis original):

La esencia de la rigidez radica en que la solución que se va a calcular varía lentamente, pero existen perturbaciones que se amortiguan rápidamente.

Shampine ^{[ 11 ]} expresa una visión similar, pero es conceptualmente más específica, relacionando la rigidez con la disipación y lo que en la práctica es la irreversibilidad :

Una forma que preferimos para describir esta última condición es que [la solución] es muy inestable en la dirección opuesta [del tiempo].

La inestabilidad inversa se puede caracterizar fácilmente utilizando normas logarítmicas . Por lo tanto, para${\displaystyle {\dot {u}}=Au}$, la norma de la solución puede acotarse utilizando desigualdades diferenciales . Sea${\displaystyle \|u\|_{2}^{2}=u^{*}u}$ser la norma euclidiana, posiblemente con${\displaystyle u\in {\mathbb {C} }^{d}}$. Entonces

${\displaystyle m_{2}[A]\cdot \|u\|_{2}\,\leq \,{\mathrm {D} }_{t}\,\|u\|_{2}\,\leq \,M_{2}[A]\cdot \|u\|_{2}\,,}$

dónde${\displaystyle {\mathrm {D} }_{t}}$denota la derivada temporal y${\displaystyle m_{2}[A]}$y${\displaystyle M_{2}[A]}$son las normas logarítmicas inferior y superior de${\displaystyle A}$, respectivamente. Estos son los extremos de una forma cuadrática simétrica,

${\displaystyle m_{2}[A]\,\leq \,{\frac {u^{*}{\mathrm {He} }(A)u}{u^{*}u}}\,\leq \,M_{2}[A]\,,}$

dónde${\displaystyle {\mathrm {He} }(A)=(A+A^{*})/2}$denota la parte hermitiana de${\displaystyle A}$Sus puntos estacionarios se obtienen a partir del problema de valores propios simétrico.

${\displaystyle {\mathrm {He} }(A)v=\mu v\,,}$

cuyo${\displaystyle d}$valores propios reales${\displaystyle \mu _{1}\geq \mu _{2}\geq \dots \geq \mu _{d}}$se denominan valores logarítmicos de${\displaystyle A}$. En relación con los valores propios${\displaystyle \lambda _{k}}$de${\displaystyle A}$, los valores logarítmicos satisfacen

${\displaystyle m_{2}[A]\,\equiv \,\mu _{d}\,\leq \,\min _{k}{\mathrm {Re} }\,\lambda _{k}\,;\qquad \max _{k}{\mathrm {Re} }\,\lambda _{k}\,\leq \,\mu _{1}\,\equiv \,M_{2}[A]\,,}$

junto con la identidad de rastreo

${\displaystyle {\mathrm {Re} }\,{\mathrm {trace} }\,A\,=\,\sum _{k}{\mathrm {Re} }\,\lambda _{k}\,=\,\sum _{k}\mu _{k}\,=\,{\mathrm {trace} }\,{\mathrm {He} }(A)\,.}$

Desde${\displaystyle u(t)={\mathrm {e} }^{tA}u_{0}}$De la desigualdad diferencial anterior se deduce que

${\displaystyle {\mathrm {e} }^{tm_{2}[A]}\,\leq \,\|{\mathrm {e} }^{tA}\|_{2}\,\leq \,{\mathrm {e} }^{tM_{2}[A]}\,\qquad t\geq 0\,,}$

limitando las tasas máximas de crecimiento y decaimiento en el sistema. Por lo tanto, durante un intervalo de tiempo de longitud${\displaystyle T>0}$, el crecimiento máximo posible en el tiempo hacia adelante es${\displaystyle {\mathrm {e} }^{TM_{2}[A]}}$, mientras que el crecimiento máximo del tiempo inverso es${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}}$.

Söderlind et al. ^{[ 12 ]} ahora cuantifican la observación de Shampine. Por lo tanto, en un sistema rígido,

${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}\,\gg \,{\mathrm {e} }^{TM_{2}[A]}\,,}$

de lo cual se deduce que${\displaystyle {\mathrm {e} }^{T{\big (}m_{2}[A]+M_{2}[A]{\big )}}\ll 1}$, es decir,${\displaystyle T{\big (}m_{2}[A]+M_{2}[A]{\big )}\ll -1}$. El indicador de rigidez se define ^{[ 13 ]}

${\displaystyle s_{2}[A]\,=\,{\frac {m_{2}[A]+M_{2}[A]}{2}}\,,}$

y una ecuación rígida se caracteriza por la condición

${\displaystyle s_{2}[A]\cdot T\ll -1\,.}$

El motivo de incluir${\displaystyle M_{2}[A]}$en la definición del indicador de rigidez es que a menudo puede ocurrir en sistemas no lineales que${\displaystyle M_{2}[A]}$es positivo, compensando la rigidez. Esto se tiene en cuenta definiendo el indicador de rigidez de manera que tenga paridad impar , es decir,${\displaystyle s_{2}[-A]=-s_{2}[A]}$.

El criterio de Shampine implica que${\displaystyle m_{2}[A]}$es grande y negativo, es decir, el límite inferior en${\displaystyle \|{\mathrm {e} }^{tA}\|_{2}}$es extremadamente pequeño. En otras palabras, el flujo es casi un semigrupo y${\displaystyle {\mathrm {e} }^{tA}}$es "casi singular". A diferencia de un sistema de Lipschitz, donde teóricamente se puede resolver la ecuación diferencial también en tiempo inverso, esto es prácticamente imposible en un sistema rígido. La ecuación rígida en tiempo inverso corresponde a un "problema inverso" extremadamente mal condicionado.

Utilizando normas logarítmicas para aplicaciones no lineales ^{[ 14 ],} los mismos argumentos se aplican a un campo vectorial no lineal.${\displaystyle f}$pero desde entonces${\displaystyle s_{2}[f]}$se define globalmente, el indicador de rigidez se define localmente como${\displaystyle s_{2}[f'(u)]}$a lo largo de la trayectoria.

Además del indicador de rigidez, Söderlind et al. ^{[ 15 ]} introducen una escala de tiempo de referencia local.${\displaystyle \Delta t(u)}$definido por

${\displaystyle \Delta t(u)\,=\,1/\max {\big (}1/T,-s_{2}[f'(u)]{\big )}\,.}$

La rigidez se caracteriza entonces por la condición independiente del método.${\displaystyle T/\Delta t(u)\gg 1}$, denominado factor de rigidez (local) . Un factor de rigidez grande es una condición necesaria para la rigidez. Pero también se pueden evaluar criterios dependientes del método. En los cálculos numéricos,${\displaystyle \Delta t(u)}$se puede comparar con el tamaño real del paso${\displaystyle h}$, cuyo factor de rigidez del tamaño del paso local${\displaystyle h/\Delta t(u)}$depende de la elección del método de integración y del requisito de precisión. Solo los métodos implícitos pueden usar tamaños de paso con${\displaystyle h/\Delta t(u)\gg 1}$, pero también puede haber partes de la integración donde el factor de rigidez del tamaño del paso sea solo moderado.

Un aspecto importante de los métodos de avance temporal implícitos es la necesidad de resolver ecuaciones de la forma

${\displaystyle u=\gamma hf(u)+\psi }$

para${\displaystyle u}$en cada paso. Aquí${\displaystyle \gamma >0}$es una constante característica del método de tamaño moderado, y${\displaystyle \psi }$es un vector conocido. Surge la pregunta de si esta ecuación tiene una solución (única) cuando${\displaystyle hf}$es "grande", en el sentido de que${\displaystyle hs_{2}[f']\ll -1}$, que representa un cálculo rígido. Reescribiendo la ecuación como

${\displaystyle {\big (}\gamma hf-I{\big )}(u)=-\psi \,,}$

Se aplica el teorema de monotonicidad uniforme: la ecuación tiene una solución única si se cumple la condición.${\displaystyle \gamma hf-I}$es monótono (negativo), es decir, si${\displaystyle M_{2}[\gamma hf-I]<0}$. Desde

${\displaystyle M_{2}[\gamma hf-I]<0\quad \Leftrightarrow \quad M_{2}[\gamma hf]<1\,,}$

Esta condición suele cumplirse con un margen razonable en cálculos rígidos.${\displaystyle M_{2}[f]<0}$, está satisfecho para todos${\displaystyle h>0}$, pero si${\displaystyle M_{2}[f]>0}$, puede requerir una restricción menor del tamaño del paso. En el criterio de Shampine, así como en el indicador de rigidez, se sostiene que

${\displaystyle {\frac {m_{2}[\gamma hf]+M_{2}[\gamma hf]}{2}}\,=\,\gamma h\!\cdot \!s_{2}[f]\ll -1\,,}$

dónde${\displaystyle m_{2}[\gamma hf]\ll -1}$y dónde${\displaystyle M_{2}[\gamma hf]}$es moderado si es positivo. Por lo tanto, aunque la ecuación normalmente requiere iteración de Newton en el cálculo rígido (${\displaystyle \gamma hf}$(no es una contracción), la existencia y la unicidad suelen estar garantizadas también cuando el factor de rigidez del tamaño del paso es grande.

Porque es relativamente costoso calcular el indicador de rigidez resolviendo el problema de valores propios simétrico.${\displaystyle \,{\mathrm {He} }{\big (}f'(u){\big )}v=\mu v\,}$para${\displaystyle \mu _{1}}$y${\displaystyle \mu _{d}}$, una alternativa económica es de interés. Observando que el indicador de rigidez es el promedio aritmético del valor logarítmico más grande y el más pequeño, una opción es reemplazar este promedio por el promedio de todos los valores logarítmicos. Con la identidad de traza mencionada anteriormente (ver arriba), definimos el funcional de traza real escalado (promedio aritmético) ^{[ 16 ].}

${\displaystyle \tau [A]\,=\,{\frac {1}{d}}\,{\mathrm {Re} }\,{\mathrm {trace} }(A)\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,a_{kk}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,\lambda _{k}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}\mu _{k}\,,}$

señalando que el cálculo de${\displaystyle \tau [A]}$es a la vez económico y trivial: no se requieren cálculos de valores propios . La traza escalada (así como la divergencia) es un funcional lineal, por lo tanto${\displaystyle \tau [-A]=-\tau [A]}$, compartiendo la extraña paridad de${\displaystyle s_{2}[A]}$. Además, para un mapa no lineal,

${\displaystyle \tau [f'(u)]\,=\,{\frac {1}{d}}\,{\mathrm {trace} }\,f'(u)\,=\,{\frac {1}{d}}\,{\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)\,=\,{\frac {1}{d}}\,({\mathrm {div} }_{u}\,f)(u)\,.}$

El funcional${\displaystyle \,\tau [f'(u)]\,}$sirve como un indicador de rigidez alternativo. Por lo tanto, la rigidez también puede caracterizarse mediante la condición de divergencia.

${\displaystyle {\frac {T}{d}}\,{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\,\ll \,-1\,,}$

lo cual se evalúa fácilmente para la mayoría de los problemas.

Para un sistema lineal${\displaystyle {\dot {u}}=Au}$, sostiene que

${\displaystyle {\mathrm {D} }_{t}\,{\mathrm {det} }\,{\mathrm {e} }^{tA}\,=\,{\mathrm {trace} }\,A\cdot {\mathrm {det} }\,{\mathrm {e} }^{tA}\,.}$

Definición de un determinante absoluto escalado (promedio geométrico)${\displaystyle \,\delta [A]=|{\mathrm {det} }\,A|^{1/d}\,}$tenemos${\displaystyle \,\delta [I]=1\,}$y${\displaystyle \,\delta [\alpha A]=|\alpha |\!\cdot \!\delta [A]\,}$, evitando la dependencia directa del determinante estándar de la dimensión${\displaystyle d}$(Esto es de especial importancia cuando la ecuación diferencial se deriva de una discretización por el método de líneas de una ecuación diferencial parcial). Aquí${\displaystyle \,\delta [A]\,}$es el "tamaño" lineal del volumen del espacio de fase abarcado por los vectores columna de${\displaystyle A}$, como en${\displaystyle \,\delta [2I]=2\,}$, independientemente de${\displaystyle d}$La ecuación diferencial para el determinante del flujo ahora se puede reescribir.

${\displaystyle {\mathrm {D} }_{t}\,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,\tau [A]\cdot \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,.}$

Por eso

${\displaystyle \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,{\mathrm {e} }^{t\,\tau [A]}\,,}$

expresando que un volumen de fase inicial (lineal) cambia por un factor${\displaystyle \,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,}$con el tiempo${\displaystyle t}$. Si${\displaystyle \,\tau [A]<0\,}$El volumen de la fase se comprime por el flujo.${\displaystyle \,{\mathrm {e} }^{tA}\,}$para${\displaystyle t\geq 0}$.

Resulta que${\displaystyle \,\tau [A]\leq 0\,}$es una condición necesaria para la estabilidad . (Comparar el criterio suficiente${\displaystyle \,M_{2}[A]\leq 0\,}$para la estabilidad de la solución cero.) La caracterización de rigidez de Shampine, en términos de soluciones fuertemente inestables en tiempo inverso, se puede escribir como${\displaystyle \,T\!\cdot \tau [-A]\gg 1\,}$, que, debido a la paridad impar de la traza, es equivalente a la condición de rigidez${\displaystyle T\!\cdot \tau [A]\ll -1}$Estas observaciones, conceptos y construcciones proporcionan la clave para caracterizar y cuantificar la rigidez, con el criterio de divergencia.${\displaystyle \,T\!\cdot \tau [f']\ll -1\,}$siendo el más simple.

Se ha recopilado una gran cantidad de problemas de prueba no triviales en el conjunto de pruebas de Bari para solucionadores de problemas de valor inicial. ^{[ 17 ]} Este contiene ejemplos detallados de problemas, solucionadores y resultados obtenidos bajo condiciones específicas por códigos bien conocidos.

Aquí se seleccionan tres problemas no lineales rígidos para ilustrar cómo se utiliza el indicador de rigidez (simplificado) para identificar y caracterizar correctamente la rigidez. Esto incluye la capacidad de distinguir el comportamiento complejo y de rápida variación en puntos de inflexión (localmente inestables), como los observados en las ecuaciones de van der Pol y Oregonator. Además, el indicador de rigidez también distingue entre subintervalos donde la rigidez es prominente y donde el problema no lo es. Finalmente, en algunos casos, incluso es posible identificar cuál de las variables dependientes contribuye más a la rigidez.

También se ilustran criterios operativos, como la dependencia del paso de tiempo adaptativo con los requisitos de precisión y su variación a lo largo de la solución, cuantificando el factor de rigidez del tamaño del paso. En conjunto, el indicador de rigidez permite obtener información detallada, incluso antes de que comience el cálculo.

Problema 1. (Propagación de la llama) Considere el modelo escalar de propagación de la llama ^{[ 18 ]} para${\displaystyle t\in [0,1]}$,

${\displaystyle \varepsilon \,{\dot {u}}\,=\,2\,(u^{2}-u^{3})\,;\qquad u(0)=\varepsilon \ll 1\,.}$

Los dos términos del lado derecho reflejan el hecho de que la combustión dentro de la llama requiere oxígeno, suministrado a través de la superficie de la llama; como${\displaystyle u}$aumenta la superficie${\displaystyle \sim u^{2}}$Con el tiempo, se vuelve demasiado pequeño para mantener la combustión en un volumen mayor.${\displaystyle \sim u^{3}}$Este modelo simplificado se utiliza habitualmente como problema de prueba para poner a prueba la adaptabilidad del paso de tiempo del software para ecuaciones rígidas, comprobando si la transición abrupta entre una región no rígida y una rígida se resuelve correctamente.

Para${\displaystyle u\in (0,1)}$, el lado derecho es positivo y la solución${\displaystyle u(t)}$es monótonamente creciente. Comenzando en${\displaystyle u(0)=\varepsilon }$, se aleja del equilibrio inestable en${\displaystyle u=0}$. La divergencia del campo vectorial (indicador de rigidez) es simplemente${\displaystyle f'(u)=2\,(2u-3u^{2})/\varepsilon }$. Desde

${\displaystyle f'(\epsilon )\approx 4\,;\qquad f'(1)=-\,{\frac {2}{\varepsilon }}\ll -1\,,}$

La solución comienza siendo no rígida, cuando${\displaystyle u\ll 1}$. En${\displaystyle t\approx 0.5}$cuando${\displaystyle u\sim 1/2}$La inestabilidad se ha vuelto catastrófica, con una transición extremadamente rápida hacia el segundo equilibrio estable.${\displaystyle u=1}$.

Cuando${\displaystyle u>2/3}$, se recupera la estabilidad y la solución entra en el régimen rígido, véase la Figura 3 , donde hemos utilizado${\displaystyle \varepsilon =10^{-3}}$La ecuación es más rígida cuanto más pequeño${\displaystyle \,\varepsilon \,}$se elige y cuanto más cerca de${\displaystyle \,u=1\,}$La solución avanza.

La demostración compara un método explícito de Runge-Kutta de tercer orden con un estimador de error de segundo orden, con un método implícito de Runge-Kutta A-estable de tercer orden, también con un estimador de error de segundo orden. El método implícito no tiene ventaja en el régimen no rígido, pero en el régimen rígido puede usar pasos mucho mayores. La diferencia entre los dos métodos aumenta para valores más pequeños de${\displaystyle \,\varepsilon }$.

Problema 2. (ecuación de van der Pol) La ecuación de van der Pol se suele presentar en un par de formas alternativas, con diferentes escalas. Es un problema no lineal cuyas soluciones se aproximan a un ciclo límite. Aquí escalamos el tiempo para escribir el sistema como

${\displaystyle {\begin{aligned}{\dot {x}}&=2\kappa \cdot y\\{\dot {y}}&=2\kappa ^{2}\cdot (1-x^{2})\,y-2\kappa \cdot x\,,\end{aligned}}}$

con condiciones iniciales${\displaystyle x(0)=2}$,${\displaystyle y(0)=0}$elegido en el ciclo límite. La escala de tiempo se ha elegido de manera que el período del ciclo límite sea igual a...${\displaystyle \mathrm {O} (1)}$casi independiente del parámetro${\displaystyle \kappa }$El problema se puede resolver entonces en${\displaystyle [0,1]}$, considerando un período completo. El problema se origina en el estudio de las oscilaciones no lineales en circuitos eléctricos.

Para valores grandes de${\displaystyle \kappa }$(véase la Figura 4 , donde${\displaystyle \kappa =200}$), el ciclo límite consta de dos ramas estables, conectadas por transiciones casi discontinuas que muestran una pérdida catastrófica de estabilidad. Debido a la identidad de la traza, el indicador de rigidez${\displaystyle s_{2}[f'(u)]}$coincide con el trazado escalado${\displaystyle \tau [f'(u)]}$, dado por

${\displaystyle s_{2}[f'(u)]\,=\,\tau [f'(u)]\,=\,\kappa ^{2}\,(1-x^{2})\,,}$

dónde${\displaystyle u}$es el vector${\displaystyle (x,y)^{\mathrm {T} }}$El indicador de rigidez se identifica fácilmente en la segunda ecuación del sistema de van der Pol, aunque pasa desapercibido.

La traza escalada tiene un término constante,${\displaystyle \kappa ^{2}}$, que representa la contribución de la parte lineal del sistema. Como este término es positivo, es desestabilizador. El equilibrio${\displaystyle x=y=0}$es inestable, ya que${\displaystyle \tau [f'(0)]=\kappa ^{2}>0}$. La traza escalada también tiene un término negativo${\displaystyle -\kappa ^{2}x^{2}}$representando la contribución no lineal a la rigidez. Como la traza escalada es independiente de${\displaystyle y}$Esta última variable no tiene ninguna influencia particular sobre la rigidez.

La solución entra en el régimen rígido cuando${\displaystyle \tau [f'(u)]\ll -1}$. Desde${\displaystyle \tau [f'(u)]=\kappa ^{2}\,(1-x^{2})}$De ello se deduce que la rigidez solo se produce para${\displaystyle \,|x|>1}$, es decir, a lo largo de las dos ramas estables, donde la rigidez es proporcional a${\displaystyle \,\kappa ^{2}}$y puede ser arbitrariamente grande. En particular, se deduce que si se utilizara un método explícito para resolver el problema a lo largo de la rama estable, el esfuerzo computacional es${\displaystyle \mathrm {O} (\kappa ^{2})}$Por el contrario, el esfuerzo es${\displaystyle \mathrm {O} (1)}$para un solucionador rígido. ^{[ 19 ]}

Cabe señalar que toda esta información está disponible a priori y se confirma mediante la solución numérica del problema.

Problema 3. (Ecuación de Oregonator) La cinética de las reacciones químicas es una fuente rica de ecuaciones diferenciales rígidas. La ecuación de Oregonator modela una reacción autocatalítica que involucra tres reacciones. Para condiciones iniciales adecuadas, tiene un ciclo límite de período${\displaystyle T\approx 305}$Las ecuaciones son

${\displaystyle {\begin{aligned}{\dot {x}}\,&=\,s\cdot (x-xy+y-q\,x^{2})\\{\dot {y}}\,&=\,(z-y-xy)/s\\{\dot {z}}\,&=\,w\cdot (x-z)\,.\end{aligned}}}$

El problema de prueba utiliza los parámetros${\displaystyle s=77.27}$,${\displaystyle q=8.375\cdot 10^{-6}}$,${\displaystyle w=0.161}$junto con las condiciones iniciales${\displaystyle x(0)=3.2,\,y(0)=1.45,\,z(0)=2.5}$. ^{[ 20 ]}

Porque la dimensión del sistema es${\displaystyle d=3}$, el indicador de rigidez${\displaystyle s_{2}[f']}$y la traza escalada (divergencia)${\displaystyle \tau [f'(u)}$no coinciden, aunque las diferencias son pequeñas. Más importante aún, dado que la ecuación es un ejemplo más complejo de rigidez, necesitamos obtener una comprensión preliminar manteniendo los parámetros el mayor tiempo posible en el análisis. Por lo tanto, elegimos la divergencia escalada, que ofrece esta ventaja. Así encontramos

${\displaystyle T\cdot \tau [f'(u)]\,=\,\alpha +\beta \,x+\gamma \,y,}$

dónde

${\displaystyle {\begin{aligned}\alpha \,&=\,T\cdot (s^{2}-sw-1)/(3s)\,\approx \,7.838\cdot 10^{3}\\\beta \,&=\,-T\cdot (2qs^{2}+1)/(3s)\,\approx \,-1.447\\\gamma \,&=\,-T\cdot s/3\,\approx \,-7.856\cdot 10^{3}.\end{aligned}}}$

Aquí el término constante${\displaystyle \alpha }$es la contribución lineal a la divergencia. Es positiva, por lo tanto desestabilizadora. Los otros dos términos representan contribuciones no lineales, ya que se multiplican${\displaystyle x}$y${\displaystyle y}$Sus coeficientes${\displaystyle \beta }$y${\displaystyle \gamma }$son negativos y contribuirán a la rigidez. Pero como la divergencia es independiente de${\displaystyle z}$, esta última variable no contribuye a la rigidez más allá del término constante.${\displaystyle \alpha }$. El coeficiente más grande es${\displaystyle \gamma }$, lo que indica que la rigidez será particularmente fuerte cuando${\displaystyle y}$es grande y positivo. Si${\displaystyle x}$Su contribución dependerá de la magnitud que alcance. En la Figura 5 , los datos reales se recopilan durante un período completo.${\displaystyle [0,305]}$graficado en el intervalo normalizado${\displaystyle [0,1]}$.

Durante un breve intervalo cuando${\displaystyle x}$es grande, la divergencia escalada es de orden${\displaystyle T\cdot \tau [f'(u)]\approx -10^{5}}$(apenas visible como una muesca en la esquina superior izquierda del gráfico de trazado escalado de alta resolución), lo que hace que el problema sea rígido allí. Pero lo peor está por venir; cuando${\displaystyle y}$a medida que aumenta la rigidez, se vuelve severa, con divergencia escalada en${\displaystyle T\cdot \tau [f'(u)]\approx -10^{7}}$.

Cuando se utiliza un método adaptativo de Runge-Kutta A-estable de tercer orden, los tamaños de paso superan la escala de tiempo de referencia en tres o cuatro órdenes de magnitud. Esto corresponde a la mejora en la eficiencia del solucionador rígido en comparación con el uso potencial de un método explícito.

Como se mencionó anteriormente, las ecuaciones rígidas requieren métodos implícitos de discretización temporal, ya sean métodos de Runge-Kutta o métodos lineales multipaso. También existen alternativas basadas en técnicas de extrapolación. Sin embargo, no todos los métodos implícitos son adecuados. Un buen método debe tener una región de estabilidad.${\displaystyle S}$cubriendo todo o la mayor parte del semiplano negativo${\displaystyle {\mathbb {C} }^{-}}$.

La región de estabilidad se determina aplicando el método a la ecuación de prueba de Dahlquist ^{[ 21 ].}${\displaystyle \,{\dot {u}}=\lambda u}$y la región de estabilidad${\displaystyle S}$es el conjunto de${\displaystyle h\lambda \in {\mathbb {C} }}$para las cuales el método produce soluciones acotadas. Un método con${\displaystyle S\supset {\mathbb {C} }^{-}}$, es decir, donde${\displaystyle S}$Contiene todo el semiplano izquierdo, se denomina A-estable.

Existen muchos métodos de Runge-Kutta A-estables de alto orden para elegir, pero desafortunadamente, para los métodos lineales multipaso A-estables el orden de convergencia está limitado a${\displaystyle p=2}$Por lo tanto, normalmente hay que conformarse con menos. Entre los métodos de pasos múltiples, la mejor opción para ecuaciones rígidas son los métodos de diferenciación hacia atrás (BDF), a los que pertenece el método implícito de Euler. Existen métodos BDF de órdenes de convergencia${\displaystyle p=1:6}$pero solo pedidos${\displaystyle p\leq 2}$son A-estables. Si bien la región de estabilidad permanece grande para órdenes superiores, eventualmente se deteriora y solo los métodos hasta el orden${\displaystyle p=5}$se utilizan con regularidad.

El software basado en métodos de tipo BDF (o similares) incluye ode15s de Matlab y códigos en C o Fortran como CVODE, LSODE, MEBDF y DASSL. Este software es complejo, robusto y de alta complejidad, utiliza adaptabilidad de orden y paso de tiempo variables, y suele ofrecer numerosas opciones para problemas con estructuras especiales o para el manejo específico de matrices jacobianas. Se pueden especificar criterios de precisión tanto para estimaciones de error absoluto como relativo. Algunos códigos incluyen opciones para ecuaciones no rígidas o para ecuaciones diferenciales-algebraicas.

Para los métodos implícitos de Runge-Kutta aplicados a la ecuación de prueba, la ecuación diferencial se reemplaza por una recursión.${\displaystyle u_{n+1}=R(h\lambda )\,u_{n}}$donde la función racional${\displaystyle R}$se denomina función de estabilidad.

${\displaystyle {\mathrm {Re} }\,z\leq 0\,\Rightarrow \,|R(z)|\leq 1\,,}$

El método es A-estable. Existen métodos A-estables de alto orden, pero pocas implementaciones. Una ventaja de los métodos A-estables de Runge-Kutta es que

${\displaystyle M_{2}[hA]\leq 0\,\Rightarrow \,\|R(hA)\|_{2}\leq 1\,,}$

es decir, si${\displaystyle A}$Si es definida negativa, entonces el método produce una recursión contractiva. (Esta es una variante de la desigualdad de von Neumann ). Desafortunadamente, esto solo se extiende a problemas no lineales bajo condiciones adicionales. Los métodos de Runge-Kutta B-estables son un subconjunto de los métodos A-estables, de modo que si${\displaystyle \,M_{2}[hf]\leq 0\,}$En un problema no lineal, el método produce una recursión contractiva. ^{[ 22 ] [ 23 ] [ 24 ]} También son comunes otros requisitos de estabilidad especiales, como la estabilidad L. ^{[ 25 ]} Este es otro subconjunto del método A-estable, con el requisito adicional${\displaystyle \,R(\infty )=0\,}$para una mejor amortiguación de los modos propios para los cuales${\displaystyle \,{\mathrm {Re} }\,h\lambda \ll -1}$.

Debido a que los métodos implícitos de Runge-Kutta tienen una alta complejidad computacional, también se consideran requisitos especiales de eficiencia cuando se desarrolla el software. Entre el software eficiente de Runge-Kutta para problemas rígidos se encuentran los ode23s de Matlab de órdenes${\displaystyle 2}$y${\displaystyle 3}$y el código Fortran RADAU5 ^{[ 26 ]} , que implementa la estabilidad B y L.${\displaystyle 5^{\mathrm {th} }}$método de orden Radau IIa. Para una revisión general de los métodos de Runge-Kutta diagonalmente implícitos utilizados para ecuaciones rígidas, véase Kennedy y Carpenter. ^{[ 27 ]}

En todos los casos, la solución numérica práctica de ecuaciones rígidas requiere software especializado y profesional. Aunque se trata de una tarea rutinaria, el éxito suele exigir un estudio minucioso de la configuración adecuada de los códigos y una comprensión precisa del problema a resolver.

1. ¿Se comprende completamente la rigidez? En la literatura, a veces se sugiere que no existe una definición precisa de rigidez. Esto parece deberse a la mención o el uso generalizado del coeficiente de rigidez, cuyas evidentes deficiencias se reconocen desde hace tiempo. En ocasiones, se recurre a criterios operacionales, como el método de discretización y los requisitos de precisión, para resolver estos problemas. Sin embargo, esto complica en exceso lo que, en esencia, es un asunto sencillo. Por lo tanto, para el profesional es obvio que existe una distinción entre ecuaciones no rígidas y rígidas; en consecuencia, debe ser posible describir matemáticamente esta distinción. El indicador de rigidez proporciona un criterio simple y necesario para caracterizar y cuantificar esta distinción. Hoy en día, se puede afirmar con seguridad que la rigidez es un fenómeno complejo, pero plenamente comprendido, y que se dispone de software especializado excelente, eficiente y fiable.

2. Indicador de rigidez. El indicador de rigidez${\displaystyle s_{2}[f']}$es más robusto que la divergencia escalada${\displaystyle \tau [f']}$Mientras que el primero solo involucra los dos valores logarítmicos extremos, el segundo también utiliza los valores logarítmicos intermedios, a pesar de que no tienen un impacto similar en la estabilidad. La ventaja de la divergencia escalada es que evita los cálculos de valores propios y permite una comprensión analítica a priori de la dinámica, como se ilustra en los ejemplos anteriores. Generalmente, una cuantificación aproximada es suficiente, ya que la rigidez suele ser una cuestión de órdenes de magnitud .

3. Sistemas disipativos frente a sistemas conservativos. La noción de rigidez, como se mencionó anteriormente, es una medida de disipatividad, es decir, de amortiguación o pérdida de energía. Algunas fuentes sugieren que también existen sistemas conservativos rígidos, en particular si se derivan de discretizaciones del método de líneas de EDP hiperbólicas . Pero estos sistemas pertenecen a una categoría diferente: la de sistemas (potencialmente) altamente oscilatorios. Con autovalores ubicados muy lejos en el eje imaginario, las altas frecuencias correspondientes pueden suprimirse eligiendo un método implícito con pasos grandes. Sin embargo, de acuerdo con el teorema de muestreo , los fenómenos de alta frecuencia ya no se resuelven y las formas de onda pueden distorsionarse. La distinción entre problemas disipativos y conservativos es significativa, y el enfoque computacional es diferente. Mientras que la rigidez está estrechamente relacionada con la irreversibilidad de los problemas parabólicos , los problemas hiperbólicos no tienen amortiguación intrínseca y presentan invariantes, como la conservación de la energía. De forma similar, los sistemas hamiltonianos separables poseen un campo vectorial libre de divergencia y, por lo tanto, conservan el área (el volumen de fase). Debido a esta estructura, el indicador de rigidez es cero.

4. Etimología. El término "rígido" fue introducido por Curtiss y Hirschfelder. Según Hirschfelder, ^{[ 28 ]} el término fue elegido porque los primeros ejemplos estaban asociados con sistemas servo, en los que existía un "acoplamiento estrecho" entre el servo y el sistema controlado. Un efecto similar también puede observarse en sistemas de control con retroalimentación negativa de alta ganancia, donde una mayor ganancia aumenta la rigidez. Otra conexión sugerida es la noción (mecánica) de rigidez, representada por la ecuación de segundo orden.

${\displaystyle m{\ddot {x}}+c{\dot {x}}+kx=F\,,}$

dónde${\displaystyle m}$representa la masa,${\displaystyle c}$el coeficiente de amortiguación,${\displaystyle k}$la constante elástica y${\displaystyle F}$una fuerza externa aplicada. La idea es que las ecuaciones con "resortes rígidos" (grandes)${\displaystyle k}$) causan rigidez en el sentido matemático. Tomando${\displaystyle m=1}$, la ecuación se reescribe como

${\displaystyle {\begin{aligned}{\dot {x}}&=y\\{\dot {y}}&=-kx-cy+F\,,\end{aligned}}}$

que es un sistema de la forma${\displaystyle {\dot {u}}=Au+g}$. Los valores propios de${\displaystyle A}$son

${\displaystyle \lambda _{1,2}\,=\,{\frac {-c\pm {\sqrt {c^{2}-4k}}}{2}}\,;}$

Se afirma que estos son grandes si${\displaystyle k}$es grande. Para amortiguamiento crítico , los parámetros deben satisfacer${\displaystyle c=2{\sqrt {mk}}}$, que con${\displaystyle m=1}$conduce a

${\displaystyle \lambda _{1,2}\,=\,-c/2\,.}$

Curiosamente, al calcular el indicador de rigidez (en este caso idéntico a la divergencia escalada) del campo vectorial se obtiene

${\displaystyle s_{2}[A]\,=\,\tau [A]\,=\,-c/2\,,}$

que es grande si y solo si la constante de amortiguación${\displaystyle c}$es grande, independientemente de la constante del resorte . Por lo tanto, la única causa de rigidez es el amortiguador, que disipa energía. Es difícil evitar la conclusión de que el término "rígido" es un nombre inapropiado, que rápidamente se estableció en un nuevo contexto con un significado diferente. Si bien en ingeniería mecánica una constante de resorte grande generalmente va acompañada de una constante de amortiguamiento grande para tener un amortiguamiento (casi) crítico, la confusión puede entenderse.

• Fórmula de diferenciación hacia atrás , una familia de métodos implícitos especialmente utilizados para la solución de ecuaciones diferenciales rígidas.
• Número de condición
• La inclusión diferencial , una extensión de la noción de ecuación diferencial que permite discontinuidades, se utiliza en parte para evitar algunos problemas de rigidez.
• Métodos explícitos e implícitos

• Introducción al modelado basado en la física: funciones de energía y rigidez.
• Sistemas rígidos Lawrence F. Shampine y Skip Thompson Scholarpedia , 2(3):2855. doi:10.4249/scholarpedia.2855