De una línea de fase a un espacio de estados

En el capítulo 33 estudiamos ecuaciones autónomas de una variable, de la forma x˙=f(x)\dot x=f(x). En ellas la recta impone restricciones muy fuertes. Entre dos equilibrios consecutivos, el signo del campo determina una dirección única de movimiento. Las trayectorias no pueden rodear un punto, describir espirales ni combinar varias escalas de evolución independientes. Al aumentar el número de variables cambia la geometría del problema.

Consideremos ahora un vector de estado x(t)=(x1(t),…,xn(t))T\mathbf x(t)=(x_1(t),\ldots,x_n(t))^{\mathsf T} y un sistema lineal autónomo con coeficientes constantes,

x˙=Ax,\boxed{\dot{\mathbf x}=A\mathbf x,}

donde A∈Rn×nA\in\mathbb R^{n\times n}. Cada componente de la derivada puede depender de varias componentes del estado. En dos dimensiones,

x˙=a11x+a12y,y˙=a21x+a22y.\dot x=a_{11}x+a_{12}y, \qquad \dot y=a_{21}x+a_{22}y.

Los términos diagonales describen cómo contribuye una variable a su propia tasa de cambio, mientras que los términos fuera de la diagonal describen los acoplamientos entre variables. Esta lectura es algebraica. Atribuir a un coeficiente una interpretación física o causal adicional exige especificar el modelo que produjo la matriz.

La palabra lineal tiene aquí un significado exacto. Si x(t)\mathbf x(t) e y(t)\mathbf y(t) son soluciones, entonces c1x(t)+c2y(t)c_1\mathbf x(t)+c_2\mathbf y(t) también es solución para constantes cualesquiera c1c_1 y c2c_2. Este principio de superposición permitirá descomponer una trayectoria en componentes dinámicas más simples, llamadas modos.

El paso desde una dimensión conserva lo aprendido en los capítulos anteriores. Siguen haciendo falta la existencia y la unicidad. También sigue siendo necesario distinguir la estabilidad de la atracción. La novedad es que la dirección del movimiento ya no puede resumirse con un signo, de modo que hay que estudiar cómo la matriz transforma las distintas direcciones del espacio de estados.

Unidades y significado de la matriz

Antes de manipular valores propios hay que comprobar qué describe realmente AA. Si todas las componentes de x\mathbf x tienen la misma unidad [X][X] y el tiempo tiene unidad [T][T], cada coeficiente de la matriz debe tener unidad [T]−1[T]^{-1}, porque [Ax]=[X]/[T][A\mathbf x]=[X]/[T]. La situación cambia si las componentes tienen unidades distintas. Si xix_i tiene unidad [Xi][X_i] y xjx_j unidad [Xj][X_j], entonces

[aij]=[Xi][Xj][T].[a_{ij}]=\frac{[X_i]}{[X_j][T]}.

Una matriz puede, por tanto, combinar coeficientes con unidades diferentes. Escribir todos sus elementos como números directamente comparables puede ocultar una inconsistencia dimensional.

Un cambio legítimo de unidades modifica a la vez el vector y la matriz. Si z=Sx\mathbf z=S\mathbf x con SS invertible, entonces z˙=SAS−1z\dot{\mathbf z}=SAS^{-1}\mathbf z. Las matrices AA y SAS−1SAS^{-1} son semejantes y tienen los mismos valores propios, aunque sus entradas y sus vectores propios escritos en coordenadas cambien. Esta invariancia será importante, porque la tasa de un modo no debería depender de que una concentración se exprese en gramos o en miligramos, siempre que el sistema se transforme de manera coherente.

No ocurre lo mismo con la apariencia geométrica de una figura. Un reescalamiento no uniforme de los ejes modifica ángulos, distancias y la forma visual de las trayectorias. Los valores propios de un operador lineal son invariantes bajo semejanza, mientras que la longitud euclídea de un vector dibujado en coordenadas arbitrarias no lo es.

Existencia global y exponencial de una matriz

El campo f(x)=Axf(\mathbf x)=A\mathbf x es globalmente Lipschitz en cualquier norma subordinada, porque

∥Ax−Ay∥≤∥A∥ ∥x−y∥.\|A\mathbf x-A\mathbf y\| \le \|A\|\,\|\mathbf x-\mathbf y\|.

Por el teorema de Picard-Lindelöf, cada condición inicial determina una solución local única. La cota ∥Ax∥≤∥A∥ ∥x∥\|A\mathbf x\|\le \|A\|\,\|\mathbf x\| da además un crecimiento lineal. Con ella, la desigualdad de Grönwall impide una explosión en tiempo finito como la de x′=x2x'=x^2 en el capítulo 32. El sistema lineal homogéneo con matriz constante tiene, por tanto, una solución global única para cualquier condición inicial finita (Hirsch et al., 2013; Teschl, 2012).

La solución general se escribe mediante la exponencial matricial,

x(t)=eA(t−t0)x0,eAt=I+At+A2t22!+A3t33!+⋯ .\boxed{\mathbf x(t)=e^{A(t-t_0)}\mathbf x_0,} \qquad e^{At}=I+At+\frac{A^2t^2}{2!}+\frac{A^3t^3}{3!}+\cdots.

La serie converge para toda matriz cuadrada y satisface

ddteAt=AeAt=eAtA,eA0=I,\frac{d}{dt}e^{At}=Ae^{At}=e^{At}A, \qquad e^{A0}=I,

de modo que la expresión anterior cumple la ecuación y la condición inicial (Higham, 2008). La familia φt(x0)=eAtx0\varphi_t(\mathbf x_0)=e^{At}\mathbf x_0 define un flujo global y verifica eA(t+s)=eAteAse^{A(t+s)}=e^{At}e^{As}.

La fórmula es exacta, aunque eso no significa que la serie deba calcularse término a término. Los algoritmos numéricos robustos para la exponencial matricial usan factorizaciones y aproximaciones diseñadas para controlar el costo y la estabilidad (Higham, 2008). En este capítulo usaremos scipy.linalg.expm como referencia computacional, junto con derivaciones analíticas independientes.

Hay una extensión útil. El sistema afín x˙=Ax+b\dot{\mathbf x}=A\mathbf x+\mathbf b no cumple el principio de superposición, pero si tiene un equilibrio x∗\mathbf x^* con Ax∗+b=0A\mathbf x^*+\mathbf b=0, el cambio u=x−x∗\mathbf u=\mathbf x-\mathbf x^* da u˙=Au\dot{\mathbf u}=A\mathbf u. Por eso el análisis de los sistemas lineales describe también las desviaciones alrededor del equilibrio de muchos modelos afines. En el capítulo siguiente servirá, además, como aproximación local de los sistemas no lineales.

Vectores propios como direcciones que no cambian de orientación

Un vector no nulo v\mathbf v es vector propio de AA si existe un escalar λ\lambda tal que Av=λvA\mathbf v=\lambda\mathbf v. Si la condición inicial está exactamente en esa dirección, con x0=cv\mathbf x_0=c\mathbf v, entonces

x(t)=ceλtv.\mathbf x(t)=ce^{\lambda t}\mathbf v.

La trayectoria permanece sobre la recta generada por v\mathbf v y solo cambia su amplitud. El valor propio λ\lambda determina la tasa temporal de ese modo. Si λ\lambda es real y negativo, la amplitud decrece exponencialmente. Si es positivo, crece. Si λ=0\lambda=0, la amplitud se mantiene constante en esa dirección, salvo que una estructura de Jordan introduzca términos adicionales. El valor propio tiene unidades de tiempo inverso, de modo que 1/∣λ∣1/|\lambda|, cuando λ≠0\lambda\neq0, es una escala temporal característica del modo.

Cuando una matriz tiene nn vectores propios linealmente independientes, puede formarse P=(v1  ⋯  vn)P=(\mathbf v_1\;\cdots\;\mathbf v_n) y escribirse A=PΛP−1A=P\Lambda P^{-1}, con Λ=diag⁡(λ1,…,λn)\Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n). Entonces

eAt=PeΛtP−1,eΛt=diag⁡(eλ1t,…,eλnt).e^{At}=Pe^{\Lambda t}P^{-1}, \qquad e^{\Lambda t} = \operatorname{diag}(e^{\lambda_1t},\ldots,e^{\lambda_nt}).

Si x0=c1v1+⋯+cnvn\mathbf x_0=c_1\mathbf v_1+\cdots+c_n\mathbf v_n, se obtiene

x(t)=∑j=1ncjeλjtvj.\boxed{ \mathbf x(t)= \sum_{j=1}^{n}c_j e^{\lambda_jt}\mathbf v_j. }

Esta es la descomposición modal. Las variables observadas son combinaciones de modos. Cada modo es, a su vez, una combinación de variables. En general no hay una correspondencia uno a uno entre una variable y una escala temporal.

Dos variables y dos velocidades de decaimiento

Consideremos

x˙=Ax,A=(−211−2),x(0)=(31),\dot{\mathbf x}=A\mathbf x, \qquad A= \begin{pmatrix} -2&1\\ 1&-2 \end{pmatrix}, \qquad \mathbf x(0)= \begin{pmatrix} 3\\1 \end{pmatrix},

es decir, x˙=−2x+y\dot x=-2x+y y y˙=x−2y\dot y=x-2y. Los valores propios cumplen

det⁡(A−λI)=(−2−λ)2−1=λ2+4λ+3=0,\det(A-\lambda I)=(-2-\lambda)^2-1=\lambda^2+4\lambda+3=0,

de modo que λ1=−1\lambda_1=-1 y λ2=−3\lambda_2=-3. Un vector propio para λ1\lambda_1 es v1=(1,1)T\mathbf v_1=(1,1)^{\mathsf T}, que describe un modo en que ambas variables cambian juntas. Para λ2\lambda_2 puede tomarse v2=(1,−1)T\mathbf v_2=(1,-1)^{\mathsf T}, que describe una diferencia entre las variables.

La condición inicial se descompone como (3,1)T=2v1+v2(3,1)^{\mathsf T}=2\mathbf v_1+\mathbf v_2, de modo que la solución es

x(t)=2e−t(11)+e−3t(1−1),\mathbf x(t)= 2e^{-t} \begin{pmatrix}1\\1\end{pmatrix} + e^{-3t} \begin{pmatrix}1\\-1\end{pmatrix},

y por componentes

x(t)=2e−t+e−3t,y(t)=2e−t−e−3t.\boxed{x(t)=2e^{-t}+e^{-3t},} \qquad \boxed{y(t)=2e^{-t}-e^{-3t}.}

La solución separa dos tiempos. El modo diferencial e−3te^{-3t} desaparece con rapidez, mientras que el modo común e−te^{-t} persiste más tiempo. La razón entre sus amplitudes es e−3t/(2e−t)=e−2t/2e^{-3t}/(2e^{-t})=e^{-2t}/2, que en t=2t=2 vale unos 0,009160{,}00916. La trayectoria ya está entonces casi alineada con v1\mathbf v_1, aunque ambas variables sigan lejos de cero (figura 34.1).

Descomposición de la evolución de un sistema simétrico en sus modos propios.
Figura 34.1. Descomposición modal del sistema simétrico. En el panel (a), las componentes x(t) e y(t) se comparan con el modo lento, proporcional a (1, 1), así como con el modo rápido, proporcional a (1, −1). En el panel (b), la trayectoria que parte de (3, 1) se alinea con la dirección lenta antes de llegar al origen. Ampliar figura ↗

El ejemplo contiene un detalle instructivo. Aunque ambos valores propios son negativos, x˙(0)=Ax0=(−5,1)T\dot{\mathbf x}(0)=A\mathbf x_0=(-5,1)^{\mathsf T}, de modo que la variable yy aumenta al comienzo. La estabilidad de los modos no obliga a que cada componente observada disminuya monótonamente, porque las variables son sumas de modos con signos diferentes.

Traza y determinante en dos dimensiones

Para una matriz de 2×22\times2 con entradas aa, bb, cc y dd, el polinomio característico es

λ2−(tr⁡A)λ+det⁡A=0,λ1,2=tr⁡A±(tr⁡A)2−4det⁡A2.\lambda^2-(\operatorname{tr}A)\lambda+\det A=0, \qquad \lambda_{1,2} =\frac{\operatorname{tr}A\pm \sqrt{(\operatorname{tr}A)^2-4\det A}}{2}.

La traza es la suma de los valores propios y el determinante, su producto, contando las multiplicidades. En el ejemplo anterior, tr⁡A=−4\operatorname{tr}A=-4 y det⁡A=3\det A=3, lo que da −1-1 y −3-3.

El discriminante Δ=(tr⁡A)2−4det⁡A\Delta=(\operatorname{tr}A)^2-4\det A indica si los valores propios de una matriz real de dimensión dos son dos reales distintos, un real repetido o un par complejo conjugado. Esa información organiza el retrato de fases, aunque no reemplaza el cálculo de los vectores propios ni describe toda la geometría transitoria. Dos matrices pueden compartir la traza y el determinante y tener vectores propios con orientaciones distintas.

Valores propios complejos, oscilación y envolvente exponencial

Una matriz real puede tener valores propios complejos. Como sus coeficientes son reales, aparecen en pares conjugados, λ1,2=α±iω\lambda_{1,2}=\alpha\pm i\omega. La parte real α\alpha controla una envolvente exponencial y la parte imaginaria ω\omega introduce la oscilación. El caso canónico

A=(α−ωωα)A= \begin{pmatrix} \alpha&-\omega\\ \omega&\alpha \end{pmatrix}

puede resolverse exactamente. Para x(0)=(1,0)T\mathbf x(0)=(1,0)^{\mathsf T},

x(t)=eαtcos⁡(ωt),y(t)=eαtsin⁡(ωt),x(t)=e^{\alpha t}\cos(\omega t), \qquad y(t)=e^{\alpha t}\sin(\omega t),

y la distancia euclídea al origen en estas coordenadas canónicas es eαte^{\alpha t}. Si α<0\alpha<0 la trayectoria describe una espiral hacia el origen, si α>0\alpha>0 se aleja de él y si α=0\alpha=0 recorre una circunferencia de radio constante. El período es T=2π/∣ω∣T=2\pi/|\omega|.

Con α=−0,2\alpha=-0{,}2 y ω=2\omega=2, los valores propios son −0,2±2i-0{,}2\pm2i. La amplitud decae con una constante temporal de 5 y el período es π\pi, de modo que en cada vuelta el radio queda multiplicado por e−0,2π≈0,533e^{-0{,}2\pi}\approx0{,}533 (figura 34.2).

Sistema canónico con valores propios −0,2 ± 2i desde (1, 0).
Figura 34.2. Sistema canónico con valores propios −0,2 ± 2i desde (1, 0). La trayectoria describe una espiral de período π y su radio disminuye como e^(−0,2t). Las circunferencias marcan el radio inicial y el radio al cabo de un período. Ampliar figura ↗

La interpretación debe matizarse para una matriz general con el mismo par complejo. Una transformación lineal puede llevar el sistema a esta forma sobre el plano invariante, pero esa transformación no tiene por qué conservar distancias ni ángulos. En las coordenadas originales pueden aparecer elipses o espirales deformadas. Lo que queda ligado al espectro es la tasa exponencial asintótica y la frecuencia del modo. Una circunferencia euclídea visible depende, en cambio, de las coordenadas.

Cuando no existe una base completa de vectores propios

La descomposición en nn modos propios independientes exige suficientes vectores propios. No toda matriz es diagonalizable. Consideremos

A=(−110−1).A= \begin{pmatrix} -1&1\\ 0&-1 \end{pmatrix}.

El único valor propio es −1-1, con multiplicidad algebraica dos, pero el espacio propio tiene dimensión uno y la matriz contiene un bloque de Jordan. Puede escribirse A=−I+NA=-I+N, con N=(0100)N=\begin{pmatrix}0&1\\0&0\end{pmatrix} y N2=0N^2=0. Como −I-I y NN conmutan,

eAt=e−teNt=e−t(I+tN)=e−t(1t01).e^{At}=e^{-t}e^{Nt} =e^{-t}(I+tN) =e^{-t} \begin{pmatrix} 1&t\\0&1 \end{pmatrix}.

Desde x0=(0,1)T\mathbf x_0=(0,1)^{\mathsf T} se obtiene x(t)=e−t(t,1)T\mathbf x(t)=e^{-t}(t,1)^{\mathsf T}. Aparece un factor polinómico tt que multiplica a la exponencial. La trayectoria sigue tendiendo a cero porque te−t→0te^{-t}\to0, aunque no puede escribirse como suma de dos modos asociados a dos vectores propios distintos.

En bloques de Jordan mayores aparecen potencias superiores de tt, lo que importa para formular correctamente la estabilidad. Con parte real negativa, los factores polinómicos no impiden el decaimiento asintótico. Si el valor propio está sobre el eje imaginario, en cambio, un bloque de Jordan no trivial puede producir un crecimiento polinómico que destruye la estabilidad de Lyapunov (Hirsch et al., 2013; Teschl, 2012).

El espectro y la estabilidad del origen

En el sistema lineal x˙=Ax\dot{\mathbf x}=A\mathbf x, el origen es siempre un equilibrio. La relación entre el espectro y su estabilidad es especialmente precisa (Hirsch et al., 2013; Perko, 2001; Strogatz, 2015). Si todos los valores propios cumplen Re⁡(λj)<0\operatorname{Re}(\lambda_j)<0, el origen es globalmente exponencialmente estable. Existe entonces una cota del tipo

∥eAt∥≤Me−γt\|e^{At}\|\le Me^{-\gamma t}

con constantes apropiadas M≥1M\ge1 y γ>0\gamma>0. Todas las trayectorias convergen entonces a cero. Si al menos un valor propio cumple Re⁡(λj)>0\operatorname{Re}(\lambda_j)>0, el origen es inestable, porque hay una dirección o un subespacio donde las perturbaciones crecen exponencialmente.

El caso frontera requiere más cuidado. Si todos los valores propios tienen parte real no positiva y los de parte real nula son semisimples, es decir, no pertenecen a bloques de Jordan no triviales, el origen puede ser estable sin ser asintóticamente estable. El ejemplo más simple es una rotación pura con valores propios ±i\pm i, cuyas trayectorias permanecen acotadas sin converger al origen. Si un valor propio del eje imaginario tiene un bloque de Jordan de tamaño mayor que uno, aparecen factores polinómicos y el origen es inestable.

Esta clasificación es exacta para los sistemas lineales de coeficientes constantes. En el capítulo 35 usaremos el espectro del jacobiano para estudiar un sistema no lineal cerca de un equilibrio. Allí las conclusiones tendrán un alcance local y aparecerán casos en que la linealización no decide.

Amplificación transitoria con valores propios negativos

Hay una simplificación tentadora según la cual, si todos los valores propios tienen parte real negativa, cualquier medida razonable de la perturbación disminuye desde el primer instante. Esa conclusión es falsa. Consideremos

A=(−180−2),x(0)=(01).A= \begin{pmatrix} -1&8\\ 0&-2 \end{pmatrix}, \qquad \mathbf x(0)= \begin{pmatrix} 0\\1 \end{pmatrix}.

Los valores propios son −1-1 y −2-2, de modo que el origen es asintóticamente estable. Los vectores propios pueden elegirse como v1=(1,0)T\mathbf v_1=(1,0)^{\mathsf T} y v2=(−8,1)T\mathbf v_2=(-8,1)^{\mathsf T}. La condición inicial contiene una cancelación grande, (0,1)T=8v1+v2(0,1)^{\mathsf T}=8\mathbf v_1+\mathbf v_2. Los dos modos decaen a velocidades distintas. Cuando el segundo desaparece más rápido que el primero, la cancelación inicial se rompe. La solución exacta es

x(t)=8(e−t−e−2t),y(t)=e−2t.\boxed{x(t)=8(e^{-t}-e^{-2t}),} \qquad \boxed{y(t)=e^{-2t}.}

Aunque ∥x(0)∥2=1\|\mathbf x(0)\|_2=1, la norma baja apenas a 0,969 cerca de t≈0,03t\approx0{,}03 y luego crece hasta unos 2,016 cerca de t≈0,677t\approx0{,}677, antes de comenzar su decaimiento definitivo (figura 34.3).

Amplificación transitoria en un sistema no normal con valores propios −1 y −2.
Figura 34.3. Amplificación transitoria en un sistema no normal con valores propios −1 y −2. La norma euclídea de la trayectoria que parte de (0, 1) supera el doble de su valor inicial antes de decaer, porque la componente x crece mientras se deshace la cancelación entre los dos modos. Ampliar figura ↗

La matriz no es normal, porque no conmuta con su traspuesta, ATA≠AATA^{\mathsf T}A\neq AA^{\mathsf T}. Sus vectores propios tampoco son ortogonales. En los sistemas no normales, los valores propios determinan el comportamiento asintótico, pero pueden describir de forma incompleta la amplificación a tiempo finito. El problema es central en el análisis de la estabilidad hidrodinámica y en la teoría de los pseudoespectros (Trefethen y Embree, 2005).

La lección es precisa. El espectro responde muy bien a la pregunta por la estabilidad asintótica lineal, aunque no agota las preguntas sobre el crecimiento transitorio. Una variable, una norma o un observable pueden aumentar temporalmente dentro de un sistema cuyo equilibrio es globalmente asintóticamente estable.

Modos, variables y observables

El ejemplo anterior muestra otra distinción importante. Un modo propio es una dirección del operador lineal, una variable de estado es una coordenada elegida para describir el sistema y un observable puede ser, además, una combinación de esas variables.

En un modelo lineal de compartimentos, por ejemplo, dos cantidades físicas x1x_1 y x2x_2 pueden producir una curva observada que es suma de varias exponenciales. Las constantes temporales de esas exponenciales corresponden a modos del sistema completo, de modo que identificar automáticamente una exponencial con un compartimento es incorrecto. Cada modo puede involucrar varios compartimentos y cada compartimento puede contener varios modos.

La interpretación depende también de la condición inicial. Si x0\mathbf x_0 carece de componente sobre un vector propio determinado, ese modo no aparece en la trayectoria aunque forme parte del sistema. Una sola ejecución puede, por tanto, ocultar modos existentes. A la inversa, una componente observada puede mostrar un máximo, un cambio de signo o un transitorio que no corresponde a un valor propio positivo. En el ejemplo simétrico, y(t)y(t) aumenta al comienzo aunque ambos modos sean contractivos, porque el movimiento de una coordenada resulta de la interferencia entre modos.

La distinción limita las inferencias mecanísticas a partir de series temporales. Ajustar una suma de exponenciales puede sugerir escalas temporales, pero identificar la estructura interna que las produjo requiere información adicional sobre la observación, las entradas, los parámetros y la identificabilidad.

Implementación y comprobación computacional

La descomposición modal puede verificarse con pocas operaciones. El bloque siguiente calcula los valores propios del ejemplo simétrico, compara la exponencial matricial con la solución cerrada y localiza el máximo de la norma en el ejemplo no normal.

import numpy as np
from scipy.linalg import expm

A = np.array([[-2.0, 1.0],
              [1.0, -2.0]])
x0 = np.array([3.0, 1.0])

valores, vectores = np.linalg.eigh(A)
print(valores)

for t in (0.0, 1.0, 2.0):
    exacta = np.array([2*np.exp(-t) + np.exp(-3*t),
                       2*np.exp(-t) - np.exp(-3*t)])
    print(t, np.max(np.abs(expm(A*t) @ x0 - exacta)))

B = np.array([[-1.0, 8.0],
              [0.0, -2.0]])
tiempos = np.linspace(0.0, 5.0, 5001)
normas = [np.linalg.norm(expm(B*s) @ np.array([0.0, 1.0])) for s in tiempos]
k = int(np.argmax(normas))
print(tiempos[k], normas[k])

La salida da primero los valores propios −3-3 y −1-1. Las diferencias entre la exponencial matricial y la solución cerrada no superan 3×10−153\times10^{-15}. La última línea localiza el máximo de la norma en t=0,677t=0{,}677, con valor 2,0162{,}016. Como AA es simétrica, eigh devuelve vectores propios ortonormales y valores propios reales. Para una matriz general habría que usar eig, sin suponer que los vectores propios sean ortogonales ni que estén bien condicionados. Estas cifras verifican los cálculos del capítulo y no aportan evidencia empírica sobre ningún sistema real.

Errores frecuentes al interpretar sistemas lineales

Un error habitual consiste en leer cada valor propio como si perteneciera a una variable. Los valores propios pertenecen al operador completo. Cada variable coincide con un modo independiente solo cuando la matriz ya es diagonal en las coordenadas elegidas. También es un error afirmar que los valores propios negativos obligan a todas las componentes a decrecer monótonamente. El ejemplo simétrico muestra una componente que crece al comienzo. El ejemplo no normal muestra incluso un crecimiento transitorio de la norma total.

La parte real y la parte imaginaria de un valor propio complejo responden a preguntas distintas. La parte real controla el crecimiento o el decaimiento exponencial, mientras que la imaginaria determina la frecuencia del modo, de modo que descartar cualquiera de ellas empobrece la descripción. Un valor propio repetido, por otra parte, no garantiza dos vectores propios, porque la multiplicidad algebraica y la dimensión del espacio propio son conceptos diferentes. Una matriz defectiva produce factores polinómicos además de exponenciales.

Los ángulos y las longitudes del plano de fases dependen de las unidades. Un reescalamiento de las variables puede cambiar drásticamente la apariencia de una trayectoria. Los valores propios sobreviven a un cambio lineal invertible de coordenadas por semejanza, mientras que las distancias euclídeas carecen en general de esa invariancia.

La estabilidad asintótica tampoco excluye una respuesta transitoria. Un sistema puede amplificar temporalmente una perturbación y converger después a su equilibrio. En preguntas de seguridad, de control o de respuesta máxima, el transitorio puede importar más que el límite cuando t→∞t\to\infty. Un buen ajuste lineal, por último, no demuestra que el fenómeno real sea lineal. Los sistemas lineales son modelos exactos en algunos contextos y aproximaciones en otros. El capítulo siguiente mostrará cómo un sistema no lineal puede aproximarse localmente mediante su jacobiano.

Consolidación

El paso de una a varias variables introduce una idea organizadora central, según la cual una trayectoria puede entenderse como superposición de modos de evolución. La solución global de x˙=Ax\dot{\mathbf x}=A\mathbf x es x(t)=eAtx0\mathbf x(t)=e^{At}\mathbf x_0. Cuando AA es diagonalizable, la condición inicial se descompone en vectores propios y cada componente modal evoluciona con un factor eλte^{\lambda t}. Los valores propios reales negativos producen decaimiento y los positivos, crecimiento. Los pares complejos introducen oscilaciones cuya envolvente depende de la parte real.

La descomposición explica por qué una misma trayectoria contiene varias escalas temporales. En el ejemplo con valores propios −1-1 y −3-3, el modo rápido desaparece y deja progresivamente visible el modo lento. Las variables observadas son combinaciones de esos modos. La diagonalización tampoco es universal, porque los bloques de Jordan generan factores polinómicos que multiplican a las exponenciales. Por eso los valores propios del eje imaginario obligan a examinar su estructura algebraica antes de concluir estabilidad.

En los sistemas lineales de coeficientes constantes, valores propios con parte real estrictamente negativa implican estabilidad exponencial global. La afirmación es asintótica. Las matrices no normales pueden producir una amplificación transitoria considerable con un espectro completamente estable. El capítulo 35 reemplazará la matriz fija del modelo lineal por el jacobiano de un sistema no lineal evaluado en un equilibrio. Los modos lineales describirán entonces la dinámica local. Habrá que precisar cuándo esa descripción permite inferir la estabilidad del sistema original.

Ejercicios

Ejercicio 34.1. Para A=(01−2−3)A=\begin{pmatrix}0&1\\-2&-3\end{pmatrix} y x(0)=(1,0)T\mathbf x(0)=(1,0)^{\mathsf T}, (a) calcule los valores propios y un vector propio para cada uno. (b) Descomponga la condición inicial en modos y escriba la solución. (c) ¿Qué componente cambia de signo en su derivada durante la evolución?

Ejercicio 34.2. Use la traza y el determinante para clasificar el origen en (a) A=(123−4)A=\begin{pmatrix}1&2\\3&-4\end{pmatrix} y (b) A=(−1−41−1)A=\begin{pmatrix}-1&-4\\1&-1\end{pmatrix}. En el caso oscilatorio, calcule el período y la constante temporal de la envolvente.

Ejercicio 34.3. Con S=diag⁡(1; 1000)S=\operatorname{diag}(1;\ 1000), (a) calcule SAS−1SAS^{-1} para la matriz simétrica del capítulo. (b) Compruebe que sus valores propios siguen siendo −1-1 y −3-3. (c) ¿Qué ocurre con el vector propio lento y con la apariencia del plano de fases?

Ejercicio 34.4. (a) Calcule eAte^{At} para A=(0100)A=\begin{pmatrix}0&1\\0&0\end{pmatrix} y escriba la solución desde (x0,v0)(x_0,v_0). (b) ¿Es estable el origen? (c) Compare el resultado con el de la matriz nula, que también tiene el valor propio doble 0.

Ejercicio 34.5. Para Aκ=(−1κ0−2)A_\kappa=\begin{pmatrix}-1&\kappa\\0&-2\end{pmatrix} y x(0)=(0,1)T\mathbf x(0)=(0,1)^{\mathsf T}, (a) muestre que x(t)=κ(e−t−e−2t)x(t)=\kappa(e^{-t}-e^{-2t}). (b) Calcule el instante y el valor del máximo de x(t)x(t). (c) ¿Cómo depende la amplificación transitoria de κ\kappa, si los valores propios no cambian?

Ejercicio 34.6. En el ejemplo simétrico, suponga que solo se observa z=x+yz=x+y. (a) Calcule z(t)z(t). (b) ¿Qué modo es invisible para ese observable? (c) ¿Qué observable permitiría detectarlo?

Consultar las soluciones de los ejercicios

Soluciones de los ejercicios

Solución 34.1. (a) El polinomio característico es λ2+3λ+2=0\lambda^2+3\lambda+2=0, con raíces −1-1 y −2-2. Vectores propios posibles son (1,−1)T(1,-1)^{\mathsf T} y (1,−2)T(1,-2)^{\mathsf T}. (b) La condición inicial es 2(1,−1)T−(1,−2)T2(1,-1)^{\mathsf T}-(1,-2)^{\mathsf T}, de modo que x(t)=2e−t−e−2tx(t)=2e^{-t}-e^{-2t} y v(t)=−2e−t+2e−2tv(t)=-2e^{-t}+2e^{-2t}. (c) La componente vv parte de cero, se vuelve negativa y regresa hacia cero. Su derivada cambia de signo en t=ln⁡2t=\ln2, donde vv alcanza su mínimo.

Solución 34.2. (a) La traza vale −3-3 y el determinante, −10-10. Un determinante negativo implica valores propios reales de signos opuestos, que son 2 y −5-5, de modo que el origen es un punto de silla inestable. (b) La traza vale −2-2, el determinante 5 y el discriminante −16-16, así que los valores propios son −1±2i-1\pm2i. El origen es estable, con trayectorias en espiral. El período es 2π/2=π2\pi/2=\pi y la constante temporal de la envolvente, 1.

Solución 34.3. (a) Se obtiene SAS−1=(−20,0011000−2)SAS^{-1}=\begin{pmatrix}-2&0{,}001\\1000&-2\end{pmatrix}. (b) Su traza sigue valiendo −4-4 y su determinante, 4−1=34-1=3, de modo que los valores propios son −1-1 y −3-3. (c) El vector propio lento pasa a ser S(1,1)T=(1,1000)TS(1,1)^{\mathsf T}=(1,1000)^{\mathsf T}, casi vertical. El plano de fases cambia drásticamente de aspecto, aunque las tasas de los modos son las mismas.

Solución 34.4. (a) Como A2=0A^2=0, eAt=I+Ate^{At}=I+At y la solución es x(t)=x0+v0tx(t)=x_0+v_0t, v(t)=v0v(t)=v_0. (b) El origen es inestable, porque una velocidad inicial arbitrariamente pequeña produce un desplazamiento que crece sin límite. (c) Con la matriz nula todo estado permanece fijo y el origen es estable. La diferencia se debe al bloque de Jordan del primer caso, que la matriz nula no tiene.

Solución 34.5. (a) Con y(t)=e−2ty(t)=e^{-2t}, la ecuación x˙=−x+κe−2t\dot x=-x+\kappa e^{-2t} con x(0)=0x(0)=0 tiene la solución κ(e−t−e−2t)\kappa(e^{-t}-e^{-2t}). (b) La derivada se anula cuando e−t=1/2e^{-t}=1/2, es decir, en t=ln⁡2≈0,693t=\ln2\approx0{,}693. El máximo vale κ/4\kappa/4. (c) La amplificación crece en proporción a κ\kappa, aunque los valores propios, −1-1 y −2-2, no dependen de κ\kappa. Con κ=8\kappa=8 el máximo de xx es 2, coherente con el máximo de la norma, de 2,016.

Solución 34.6. (a) Se obtiene z(t)=4e−tz(t)=4e^{-t}. (b) El modo rápido, proporcional a (1,−1)T(1,-1)^{\mathsf T}, suma cero en ese observable y no aparece. (c) La diferencia x−y=2e−3tx-y=2e^{-3t} lo muestra directamente. El ejemplo ilustra por qué la identificabilidad de los modos depende de qué combinación de variables se observa.

Referencias

Higham, N. J. (2008). Functions of matrices: Theory and computation. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898717778

Hirsch, M. W., Smale, S. y Devaney, R. L. (2013). Differential equations, dynamical systems, and an introduction to chaos (3.ª ed.). Academic Press.

Perko, L. (2001). Differential equations and dynamical systems (3.ª ed.). Springer. https://doi.org/10.1007/978-1-4613-0003-8

Strogatz, S. H. (2015). Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering (2.ª ed.). Westview Press.

Teschl, G. (2012). Ordinary differential equations and dynamical systems. American Mathematical Society.

Trefethen, L. N. y Embree, M. (2005). Spectra and pseudospectra: The behavior of nonnormal matrices and operators. Princeton University Press.