Equilibrium Analysis for Systems

So far, equilibrium analysis has focused on one-dimensional equations. If we introduce more equations to the system, making the system 2D, then we allow more diversity of system behaviour.

Imagine a system, represented with the following differential equations:

\[ \begin{align} \dot{x} (= \frac{dx}{dt}) &= ax+by \\ \dot{y} &= cx + dy\\ \end{align} \]

Or with the following vectors:

\[ \begin{align} \vec{x} = \begin{pmatrix} x\\ y \end{pmatrix}, \quad \dot{\vec{x}} = A \vec{x}, \quad A = \begin{pmatrix} a & b\\ c&d \end{pmatrix} \end{align} \]

Again, we define the notion of $x^*=0$ as a fixed point, with solutions visualised as trajectories on the $(x,y)$-plane, also known as the phase plane.

Example: Harmonic Oscillator

Take a spring with a mass hanging from some fixed point. If we rotate this so the fixed point is on the left hand side, and the equilibrium position of the spring is at $x=0$, then we can use the equation for a spring $F = -kx$ and derive a differential for the force:

\[ \begin{align} F &= ma \\ &= m \frac{d^2}{dt^2} x \\ \text{substitute into spring equation:}\\ m\frac{d^2}{dt^2} x &= -kx\\ m\frac{d^2}{dt^2} x -kx &= 0\\ \end{align} \]

For convenience, we can now take this equation and define some convenience notation:

\[ \omega^2 = \frac{k}{m} \quad v = \frac{dx}{dt} \implies \dot{x} = v \]

We can then define $\dot{v}$ in terms of $\omega$:

\[ \begin{align} m \frac{d^2x}{dt^2} + kx &= 0\\ m \dot{x} \frac{d}{dt} + kx &= 0\\ mv \frac{d}{dt} &= -kx\\ \frac{dv}{dt}m &= -kx\\ \dot{v}m &= -kx\\ \dot{v} &= -\frac{k}{m}x\\ &= -\omega^2 x \end{align} \]

On a plot of $x$ against $v$, we see closed orbits, or oscillations around the origin. The only fixed point of this system is at the origin, otherwise all the points will continually circle at the same distance from the origin.

We could solve this system to establish this, but we can see that solutions are on ellipses. TODO why

Therefore, for a spring with a hanging mass, with no friction or damping through other means occurring, the phase space of velocity on the $x$-axis against the position on the $y$-axis gives us the orbit of the spring as a circle:

Example: Decoupled Equations

Take another example:

\[ \begin{align} \dot{\vec{x}} &= a\vec{x}\\ A &= \begin{pmatrix} a & 0 \\ 0 & -1\\ \end{pmatrix}\\ \dot{x} &= ax\\ \dot{y} &= -y \end{align} \]

We can decouple these into a set of Malthusian growth equations, with the solutions:

\[ \begin{align} x(t) &= x_0 \exp(at)\\ y(t) &= y_0 \exp(-t) \end{align} \]

Depending on $a$, these phase portraits can differ. The direction for $y$ will always be stable, but the $x$ direction may not be. We can think of these as manifolds, with directions which attract or repel flow into an equilibrium point. We can show a stable manifold as:

\[ \{ x_0 : x(t) \to x^{FP} \text{ for } t \to \infty \} \]

and the unstable manifold as:

\[ \{ x_0 : x(t) \to x^{FP} \text{ for } t \to -\infty \} \]

This is because we are considering it attracted to the fixed point in the reverse time bounds.

Examples of Stable Nodes

This occurs when $a < 0$, otherwise they are unstable, as for a perturbation in $x$, the node will repel. We can visualise this for $a<-1, a=-1, -1<a<0$:

Non-Isolated Fixed Points

This would occur when $a=0$. For every point on the $x$ line, we'd have a fixed point at $y=0$, and there would be infinitely many of these, with any chosen $y$ converging to $y=0$ along the $x$-axis.

Saddle Points

These are the points where we have one stable and one unstable direction. In this instance, we'd have $a > 0$ which is stable in $y$, but not in $x$, so unless $x=0$, we diverge to infinity.

Stability Language

A fixed point (FP) x is

An attracting FP: all trajectories starting near $x$ will eventually approach it asymptotically
A globally attracting FP: if all trajectories will approach it asymptotically
A Liapunov stable FP: if all trajectories starting close to $x$ will remain close to $x$ for all times
A neutrally stable FP: If $x$ is Liapunov stable but not attracting.
A stable FP: $x$ is Liapunov stable and attracting
An unstable FP: $x$ is not attracting

Notice here that a point may not be Liapunov stable but might still be attracting. If we have the following flow on a circle:

\[ \dot{\theta} = 1 - \cos \theta, \theta \in [0,2\pi] \]

Then we can classify this as a half stable fixed point, as all trajectories end up in the fixed point, but may initially appear to be moving away.

Classification of Linear Systems

So far, we have seen that the $x,y$-directions show the asymptotic directions of trajectories, and straight line trajectories where there is an initial condition that starts on them and remains on them forever.

When looking for straight line solutions, we want to look for a vector with some direction and some solutions with exponential contractions and expansions.

We want to decouple the system into convenient directions such that we can understand the dynamics along these convenient directions. Given the following ansatz, we want to find the following parameters:

\[ \vec{x} (t)= \exp(\lambda t) \vec{v} | \vec{v} \ne 0 \]

Where $\vec{v}$ is the straight line direction, and $\lambda$ which is the rate of contraction. To do this, we can substitute into the differential equation:

\[ \dot{\vec{x}} = a \vec{x} \]

which is the vectorized form. If then differentiating with regards to time, we get the following:

\[ \lambda \exp(\lambda t)\vec{v} = \exp( \lambda t) A \vec{v} \]

We can then remove the $\exp(\lambda t)$ terms from both sides, yielding:

\[ \lambda \vec{v} = A \vec{v} \]

From here, we can deduce that if a straight-line solution exists then $\vec{v}$ is an eigenvector for $A$ and $\lambda$ is the eigenvalue. We can find the eigenvalue from the solutions of the characteristic equation, which is:

\[ \det(A - \lambda I) = 0 \]

This gives rise to a general solution in 2D, where $\Delta = \det(A)$, and $\tau = \text{Tr}(A)$:

\[ \lambda_{1/2} = \frac{\tau \pm \sqrt{\tau^2 - 4 \Delta}}{2} \]

Example

Given $\dot{x} = x+y$ and $\dot{y} = 4x - 2y$, with some initial conditions $(x_0, y_0) = (2,-3)$, we can substitute into $\dot{\vec{x}} = a \vec{x}$, with the following:

\[ \frac{d}{dt}\begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 4 & -2 \end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix} \]

We then try to find the characteristic equation, which occurs when we have $\det(A - \lambda I)$.

We therefore figure out the determinant, which can be done with $\det(A) = ad -bc$.

\[ \begin{align} A - \lambda I &= \begin{pmatrix} 1 & 1 \\ 4 & -2 \end{pmatrix} - \lambda \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\\ &= \begin{pmatrix} 1 & 1 \\ 4 & -2 \end{pmatrix} - \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix}\\ &= \begin{pmatrix} 1-\lambda & 1 \\ 4 & -2-\lambda \end{pmatrix}\\ \\ \det(A-\lambda I) &= ad-bc\\ &= (1-\lambda)(-2-\lambda) - 4(1)\\ &= -2 + 2\lambda - \lambda + \lambda^2 - 4\\ &= \lambda^2 + \lambda -6 \\ \lambda_{1/2} &= -1/2 \pm \sqrt{\frac{25}{4}}\\ &= 2,-3 \end{align} \]

Now we've got our values for $\lambda$, i.e., our eigenvalues, we can nnow substitute back into our $A \lambda I$ equation and multiply this through by $\begin{pmatrix}x\\y\end{pmatrix}$, setting this equal to zero:

\[ \begin{align} \begin{pmatrix} 1 - \lambda_{1/2} & 1 \\ 4 & -2 \lambda_{1/2} \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} &= \begin{pmatrix} 0\\ 0 \end{pmatrix}\\ \\ \begin{pmatrix} -1 & 1 \\ 4 & -4 \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} &= \begin{pmatrix} 0\\ 0 \end{pmatrix}\\ \begin{pmatrix} 4 & 1 \\ 4 & 1 \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} &= \begin{pmatrix} 0\\ 0 \end{pmatrix} \end{align} \]

If we solve these, we yield $v_1 = \begin{pmatrix}1\\1\end{pmatrix}$, and $v_2 = \begin{pmatrix}1\\-4\end{pmatrix}$, giving the eigenvectors for the solution.

This was with the predefined initial conditions, but as mentioned earlier, we can give a general solution with:

\[ \vec{x}(t) = c_1 \begin{pmatrix}1\\1\end{pmatrix} \exp(2t) + c_2 \begin{pmatrix}1\\-4\end{pmatrix} \exp(-3t) \]

The parameters $c_{1/2}$ need to be determined based on the initial conditions of the system, in this instance:

\[ \begin{pmatrix} 2\\ -3 \end{pmatrix} = c_1 \begin{pmatrix} 1\\ 1 \end{pmatrix} + c_2 \begin{pmatrix} 1\\ -4 \end{pmatrix} \]

We can then solve this for this set of initial conditions, which yields $c_{1/2} = 1,1$, so we have a specific solution for the equation.

Phase Portraits with Eigenvectors

Historically, we have considered those systems where each solution is somewhere about the $x$- or $y$-axis. When this is not the case, we make use of the phase portrait and plot the vectors. The crossing point is the fixed point, which may be stable, or unstable. For the example given a second ago, we would have the following phase portrait, with an unstable fixed point.

The directions of the arrows for the vectors passing through the point are given by the eigenvalues $\lambda_{1/2}$, with an eigenvalue of $<0$ being stable and $>1$ as unstable.

We can furthermore categorise the attraction properties without using the phase portrait. We can look at the eigenvalues for each eigenvector. The closer to zero the magnitude is, the slower it will converge or diverge from the fixed point. This means that for $\lambda_2 < \lambda_1 < 0$, the slow direction of convergence would be $\lambda_1$, and the fast direction would be $\lambda_2$.

Special Case: $\lambda_2 = \lambda_1$

If $\lambda = 0$, then the entire plane consists of fixed points. If we have two independent eigenvalues, then every vector is an eigenvalue, and thus we have a star node.

If we have only one eigenvalue $\ne 0$, then we have what is called a degenerate node.

Handling Complex Eigenvalues

Eigenvalues might be complex, i.e., of the form $\lambda_{1/2} = \alpha \pm i \omega$. We know that $x,y$ are linear combinations of the $c_{1/2} \exp(\alpha t) ( \cos(\omega t) + i \sin(\omega t))$, making use of Euler's formula to decompose a complex conjugate into its constituent parts.

The stability is based on the value of $\alpha$ in this intance, with $\alpha < 0$ being a stable spiral, $\alpha = 0$ being a 'center', and $\alpha > 0$ being an unstable spiral:

(Finally) Solving Harmonic Oscillator

If you now cast your mind back to the start of this page, you'll remember that we introduced the notion of the harmonic oscillator, with the following equation and initial conditions:

\[ \frac{md^2}{dt^2} x + kx = 0 \]

with $x(0) = 0$ and $\frac{dx}{dt}(0) = v_0$.

We defined the convenience functions:

$$ \begin{align} \dot{x} &= v\ \dot{v} &= - \omega^2x \end{align}

Where we have $\omega^2 = k/m$. From here, we can create the matrix $A$:

\[ A = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \]

If we then solve this for the determinant composed with the lambda ($\det(A - \lambda I) = 0$), then we get the following equation:

\[ \begin{align} \lambda^2 + \omega^2 &= 0\\ &= \pm \sqrt{-1} \omega\\ &= \pm i \omega \end{align} \]

We can then deduce the eigenvectors by substitution into $A - \lambda I$, giving:

\[ \begin{align} \begin{pmatrix} \mp i \omega & 1 \\ - \omega^2 & \mp i \omega \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} &= 0 &= \begin{pmatrix} \mp i/\omega \\ 1 \end{pmatrix} \end{align} \]

We can then solve for $x(t)$ and $v(t)$, by substitution into the equation shown earlier derived from Euler:

\[ \begin{pmatrix} x(t) \\ v(t) \end{pmatrix} = c_1 \begin{pmatrix} -\frac{i}{\omega} \\ 1 \end{pmatrix} \exp(i \omega t) + c_2 \begin{pmatrix} \frac{i}{\omega} \\ 1 \end{pmatrix} \exp(-i \omega t) \]

With initial conditions

\[ \begin{pmatrix} 0 \\ V_0 \end{pmatrix} = c_1 \begin{pmatrix} -\frac{i}{\omega} \\ 1 \end{pmatrix} + c_2 \begin{pmatrix} \frac{i}{\omega} \\ 1 \end{pmatrix} \quad \rightarrow \quad c_1 = c_2 = \frac{V_0}{2} \]

Then deriving the overall functions for $x(t)$ and $v(t)$:

\[ \begin{align} \begin{pmatrix} x(t) \\ v(t) \end{pmatrix} &= \frac{V_0}{2} \begin{pmatrix} -\frac{i}{\omega} \\ 1 \end{pmatrix} \exp(i \omega t) + \frac{V_0}{2} \begin{pmatrix} \frac{i}{\omega} \\ 1 \end{pmatrix} \exp(-i \omega t)\\ &= \begin{pmatrix} \frac{V_0}{\omega} \sin(\omega t) \\ V_0 \cos(\omega t) \end{pmatrix} \end{align} \]

We then see that we get periodic motion along $(\omega / v_0)^2 x^2 + 1/v_0^2 y^2 = 1$.

Key takeaways from this involve that the trajectory is dependent on the structure of the system and the initial conditions of the system. This shows that the trajectories of the harmonic osciallator are a family of ellipses, where the perturbation of one trajectory doesn't return us to the same trajectory.

In general, then, we have the following graph which allows us to categorise the stability of the system: