Press "Enter" to skip to content

Month: July 2023

Linear ODEs instabilities

Photo of Salvador Dalí (1904 - 1989)
Salvador Domingo Felipe Jacinto Dalí i Domènech, Marquess of Dalí of Púbol (1904 – 1989) Probably nonlinear, unstable, and nonmonotonic

This tiny post is devoted to a counter intuitive elementary phenomenon : for a linear ODE, with a time dependent matrix, stability is not always controlled by the spectrum of the matrix!

Autonomous case. It is well known that the long time behavior of a linear differential equation $$x'(t)=Ax(t)$$ in $\mathbb{R}^n$, when $A$ is diagonalizable, is dictated by the real part of the eigenvalues of $A$. Indeed, if $A=PDP^{-1}$ with $P$ complex invertible and $D$ complex diagonal, and $y(t):=P^{-1}x(t)$, then $$y'(t)=Dy(t)\quad\text{and}\quad\frac{\mathrm{d}}{\mathrm{d}t}\|y(t)\|^2=2\Re\langle Dy(t),y(t)\rangle=2\langle (\Re D)y(t),y(t)\rangle.$$ Hence the equation is (here exponentially) stable if all the eigenvalues of $A$ in $\mathbb{C}$ have a negative real part, in the sense that $x(t)\to0$ as $t\to\infty$ (exponentially fast) for all initial conditions $x(0)$. Another way to state it, using a matrix exponential, is $$x(t)=\mathrm{e}^{tA}x(0)=P\mathrm{e}^{tD}P^{-1}x(0)=P\mathrm{diag}(\mathrm{e}^{\mathrm{i}t\Im \lambda_k})\mathrm{diag}(\mathrm{e}^{t\Re\lambda_k})P^{-1}x(0).$$

If $A$ is not diagonalizable, we may use the Jordan factorization $A=PTP^{-1}$ where $T$ is upper triangular block diagonal with constant diagonal on each block. From the Jordan factorization we get the subtle decomposition $T=D+N$ where $D$ is diagonal and carries the eigenvalues of $T$, and $N$ is nilpotent and commutes with $D$, hence the exponential factorization $\mathrm{e}^{tT}=\mathrm{e}^{tD}\mathrm{e}^{tN}$ and $\mathrm{e}^{tN}$ is a polynomial in $t$. Thus the asymptotic behavior in long time of $$x(t)=\mathrm{e}^{tA}x(0)=P\mathrm{e}^{tT}P^{-1}x(0)=P\mathrm{e}^{tD}\mathrm{e}^{tN}P^{-1}x(0)$$ is dictated by the real part of the eigenvalues of $A$ : stability occurs when all the eigenvalues have a negative real part. Even if the eigenvalues of $A$ are all real and negative, then $$\frac{\mathrm{d}}{\mathrm{d}t}\|y(t)\|^2=2\Re\langle Ty(t),y(t)\rangle$$ can take positive values, see the two dimensional example below. But stability is about asymptotic analysis in long time which is not the monotonicity in time. The lack of monotonicity comes from the polynomial term $\mathrm{e}^{tN}$, related to the geometry of the eigenspaces, a non-spectral information ($n\geq2$). Of course $A$ is diagonalizable if and only if $N=0$.

If we give up block sparsity for $T$, we can get $P$ unitary, and $A=P(D+N)P^*$ is then the Schur unitary factorization, however, in this case, the nilpotent part $N$ no longer commute with the diagonal part $D$ in general, making more complicated the computation of the exponential.

The term “stability” comes from the fact that if $x(0)=0$ then $x(t)=0$ for all $t$, in other words $0$ is a stationary point, and we say that it is stable if there exists a neighborhood of $0$ such that the solution started from anywhere in this neighborhood converges in long time to $0$.

The equation $x'(t)=Ax(t)$ is called autonomous because $A$ does not depend on $t$.

Two dimensional case. The lack of monotonicity phenomenon is of geometric nature and requires at least two dimensions. Let us seek for a $2\times 2$ real matrix $T$ with all eigenvalues real and negative and for which $\langle Tx,x\rangle>0$ for some $x\in\mathbb{R}^2$. This matrix cannot be symmetric. Let us denote $\lambda_1<\lambda_2\leq0$ the eigenvalues, that we suppose distinct, and $\theta\in(0,\pi)$ the angle between the eigenvectors. Let us take $T$ upper triangular, hence the name, of the form $$T:=\begin{pmatrix} \lambda_1 & (\lambda_2-\lambda_1)\cot(\theta)\\0&\lambda_2\end{pmatrix}.$$ For all $x\in\mathbb{R}^2$, the quantity $$\langle Tx,x\rangle=\lambda_1x_1^2+\lambda_2x_2^2+(\lambda_2-\lambda_1)x_1x_2\cot(\theta)$$ reaches its maximum $$\frac{\lambda_1+\lambda_2+(\lambda_2-\lambda_1)\cot(\theta)}{2}$$ for $x=(\sin(a),\cos(a))$ with $a=(2\theta+\pi)/4$. Also, with $\rho:=\lambda_2/\lambda_1\in[0,1)$, we have $$\langle Tx,x\rangle>0\text{ for some $x$}\quad\text{if and only if}\quad\sin(\theta)<\frac{1-\rho}{1+\rho}.$$ This is impossible when $\theta=\pi/2$ ($T$ is symmetric). There exists a similar analysis for the case where $T$ has a double eigenvalue or two complex conjugate eigenvalues.

Nonautonomous case. Let us consider now the nonautonomous linear differential equation $$x'(t)=A(t)x(t),\quad\text{for which}\quad \frac{\mathrm{d}}{\mathrm{d}t}\|x(t)\|^2=2\langle A(t)x(t),x(t)\rangle.$$ If $-A(t)$ is symmetric positive definite for all $t$ then the system is stable and monotonic. If $A(t)$ is diagonalizable for all $t$ and $P(t)=P$ does not depend on $t$, then the matrices $\{A(t):t\}$ commute and the solution can be expressed using a simple exponential as $$x(t)=\mathrm{e}^{\int_0^tA(s)\mathrm{d}s}x(0)=P\mathrm{e}^{\int_0^t D(s)\mathrm{d}s}P^{-1}x(0),$$ showing that the knowledge of the spectrum of $A(t)$ for all $t$ controls the stability.

However, in general, when $P(t)$ depends on $t$, it turns out that we may have $\|x(t)\|\to+\infty$ as $t\to\infty$ while the eigenvalues of $A(t)$ are all real, negative, and constant in time! This means that the knowledge of the spectrum of $A(t)$ for all $t$ is not enough to describe the stability of the solution. In other words, the stability depends not only on the spectrum but also on the eigenspaces geometry dynamics through the matrix $P(t)$ in the Jordan or Schur factorization $$A(t)=P(t)T(t)P(t)^{-1}=P(t)(D(t)+N(t))P(t)^{-1}.$$ If we set $y(t):=P(t)^{-1}x(t)$ then $$y'(t)=P(t)^{-1}x'(t)+(P(t)^{-1})’x(t)=T(t)y(t)+(P(t)^{-1})’P(t)y(t).$$ The last term $(P(t)^{-1})’P(t)y(t)$, due to the dependency of $P(t)$ on $t$, is responsible of the instability phenomenon. This was apparently already known by Henri Poincaré (1854 – 1912) and Aleksandr Mikhailovich Lyapunov (1857 – 1918), who are at the origin of the extensive study of the stability of ODEs. The phenomenon is of geometric nature and requires a dimension $n\geq2$. Here is a concrete, simple, two dimensional counter example, due to Robert È. Vinograd : $$
A(t)=\begin{pmatrix}
6 \sin (12 t)-\frac{1}{2} 9 \cos (12 t)-\frac{11}{2} & \frac{9}{2} \sin (12 t)+6 \cos (12 t)+6 \\
\frac{9}{2} \sin (12 t)+6 \cos (12 t)-6 & -6 \sin (12 t)+\frac{9}{2} \cos (12 t)-\frac{11}{2} \\
\end{pmatrix},$$ for which the eigenvalues are $-1$ and $-10$ for all $t$, while $$x(t)=\frac{\mathrm{e}^{-13t}}{5}\begin{pmatrix}
\left(\mathrm{e}^{15 t}+4\right) \cos (6 t)-2 \left(\mathrm{e}^{15 t}-1\right) \sin (6 t) & 2 \left(\mathrm{e}^{15 t}-1\right) \cos (6 t)-\left(4 \mathrm{e}^{15 t}+1\right) \sin (6 t) \\
\left(\mathrm{e}^{15 t}+4\right) \sin (6 t)+2 \left(\mathrm{e}^{15 t}-1\right) \cos (6 t) & 2 \left(\mathrm{e}^{15 t}-1\right) \sin (6 t)+\left(4 \mathrm{e}^{15 t}+1\right) \cos (6 t)\end{pmatrix}x(0),$$ which may blow-up in long time for a well chosen initial condition $x(0)$ arbitrary close to $0$.

In order to get a generic example of $A(t)$ in dimension $n=2$, we can conjugate say an upper triangular matrix $T$ as above by a time dependent orthogonal matrix, for instance the rotation $$P(t)=\begin{pmatrix}\cos(\omega t) & -\sin(\omega t)\\\sin(\omega t) & \cos(\omega t)\end{pmatrix},\quad\omega\in\mathbb{R},$$ which gives $A(t)=P(t)TP(t)^{-1}$, in particular $T=A(0)$. We have $$P(t)=\mathrm{e}^{tQ}\quad\text{where}\quad Q=\begin{pmatrix}0 & -\omega\\\omega&0\end{pmatrix}.$$ Also $x'(t)=A(t)x(t)=P(t)TP(t)^{-1}x(t)$, and with $y(t):=P(t)^{-1}x(t)$, we get $$y'(t)=Ty(t)+(P(t)^{-1})’x(t)=(T-Q)y(t),$$ hence $y(t)=\mathrm{e}^{t(T-Q)}y(0)$, and $$x(t)=P(t)\mathrm{e}^{t(T-Q)}P(t)^{-1}x(0).$$ The contribution of $-Q$ in $\mathrm{e}^{t(T-Q)}$ is crucial for a possible blow-up. The problem is reduced to an autonomous ODE with matrix $T-Q$, which may have a positive eigenvalue even if $T$ does not. The stability is dictated by the real part of the eigenvalues of $T-Q$.

The Vinograd counter example is obtained from the two dimensional example with $\omega=-6$ in $Q$, and $\lambda_1=-10$, $\lambda_2=-1$, $\cot(\theta)=-4/3$ in $T$. The eigenvalues of $T-Q$ are $-13$ and $2$.

We could ask if the same phenomenon remains with a diagonal $T$, but the answer is no : indeed, if $T=D=\begin{pmatrix}a&0\\0 &a+b\end{pmatrix}$, then $T-Q$ has eigenvalues $\frac{1}{2}(2 a\pm\sqrt{b^2-4 \omega ^2}+b)$, the largest one $\frac{1}{2}(2 a+\sqrt{b^2-4 \omega ^2}+b)$ cannot be positive if both $a$ and $a+b$ are negative.

Conclusion. For the linear ODE $x'(t)=A(t)x(t)$, with $A(t)=P(t)(D(t)+N(t))P(t)^{-1}$, monotonicity as well as stability are not always dictated by the sole spectrum $D(t)$ : $N(t)$ acts as a perturbation of monotonicity while the dynamics of $P(t)$ acts as a perturbation of stability.

In the same spirit. It can be shown that $x'(t)=A(t)x(t)$ can be stable while $A(t)$ has for all $t$ real positive eigenvalues! We refer to the review article by Krešimir Josić and Robert Rosenbaum for more. The phenomenon says that the stability of the nonautonomous resolvent matrix equation $R'(t)=A(t)R(t)$ is not controlled solely by the spectrum of $A(t)$ for all $t$ in general. The same instability phenomenon takes place in discrete time, and corresponds to take long noncommutative words written with two matrix letters : the word may blow-up even if the letters are contractive, and this can be seen as a sort of geometric or directional resonance phenomenon. The phenomenas are related to the following topics between nonlinear aspects of linear algebra, analysis, ergodic theory, geometry, probability, and system control :

  • the difference between spectral radius and operator norm
  • the Horn problem on the spectrum of the sum of two matrices with prescribed spectrum
  • the pseudospectrum of nondiagonalizable matrices
  • the long time behavior of switched ergodic Markov dynamics
  • the May stability in mathematical biology and Girko random matrices
  • the multiplicative ergodic theorems and random walks on groups

Note. The instability phenomenon of nonautonomous linear ODEs is rarely taught in basic courses on the analysis of ODEs. Also, I have discovered that certain renowned experts in the analysis of PDEs are unaware of it and do not believe it at first glance. Yet, the phenomenon is not a surprise for people versed in the long time behavior of Markov switched dynamics, like me.

Further reading.

Related basic formulas, excellent mental exercises for your transportation time.$$\mathrm{e}^{\begin{pmatrix} a & b\\0 & a\end{pmatrix}}=\mathrm{e}^a\begin{pmatrix} 1 & b\\0 & 1\end{pmatrix}\quad\text{and}\quad\left\|\begin{pmatrix}a & b\\0 & a\end{pmatrix}\right\|_{\mathrm{op}}^2=a^2+\frac{b^2}{2}+\frac{|b|}{2}\sqrt{4a^2+b^2}.$$ Apart when $A$ is already in Jordan form, and beyond special cases, there is no explicit and general formula for $\mathrm{e}^A$ in terms of the entries of $A$. The matrices $P,D,N$ in the Jordan factorization $A=P(D+N)P^{-1}$ are in general not explicit in terms of the entries of $A$.

Photo of Aleksandr Mikhailovich Lyapunov (1857 - 1918)
Aleksandr Mikhailovich Lyapunov (1857 – 1918) maybe disappointed after the discovery of the instability phenomenon. If you have a smiling version of this photo, even a sarcastic one, please let me know!
Leave a Comment
Syntax · Style · .