Today we consider systems with constant coefficient matrices $\mathbf{A}$ whose eigenvalues are not real-valued. Being roots of the quadratic equation $\det(\mathbf{A} - \lambda\mathbf{I})=0$, these complex eigenvalues $\lambda_1, \lambda_2$ necessarily come in conjugate pairs: $\lambda_1=\lambda_2^*$. (For a given complex number $z=a+i b$, where $a$ and $b$ are real numbers, the complex conjugate of $z$ is the complex number $z^*=a-i b$.) We will especially focus on constructing real-valued solutions of systems of ODEs.
Consider
$$ \mathbf{x}' = \underbrace{\begin{bmatrix}3 & -2 \\ 4 & -1 \end{bmatrix}}_{\text{$\mathbf{A}$}}\mathbf{x}. $$The characteristic equation for $\mathbf{A}$ is
$$ \det(\mathbf{A} - \lambda\mathbf{I})=\lambda^2-2\lambda+5 = 0, $$which have the roots $\lambda_1 = 1+2i$ and $\lambda_2=1-2i$. These are the eigenvalues of $\mathbf{A}$.
We computed the eigenvector for $\lambda_1=1+2i$ to be
$$ \mathbf{u}_1 =\begin{bmatrix}1&1-i\end{bmatrix}, $$and since the matrix $\mathbf{A}$ is real-valued, the eigenvector $\mathbf{u}_2$ for $\lambda_2$ must be the complex conjugate of $\mathbf{u}_1$:
$$ \mathbf{u}_2 =\mathbf{u}_1^* = \begin{bmatrix}1&1+i\end{bmatrix}. $$Since $\mathbf{u}_1$ and $\mathbf{u}_2$ are associated with distinct eigenvalues, they are independent. This in fact implies that the solutions
$$ \mathbf{x}_{1}(t) = e^{\lambda_1 t}\mathbf{u}_1 = e^{(1+2i)t}\begin{bmatrix}1\\1-i\end{bmatrix} $$and
$$ \mathbf{x}_{2}(t) = e^{\lambda_2 t}\mathbf{u}_2 = e^{(1-2i)t}\begin{bmatrix}1\\1+i\end{bmatrix} $$form an independent set of solutions. Therefore we can build the general solution of the system of ODEs by taking linear combination of these solutions:
$$ \mathbf{x}(t)= c_1 e^{(1+2i)t}\begin{bmatrix}1&1-i\end{bmatrix}+c_2 e^{\lambda_2 t}\mathbf{u}_2 = e^{(1-2i)t}\begin{bmatrix}1&1+i\end{bmatrix}. $$Observe that $\mathbf{x}_{1}(t)=\mathbf{x}_{2}^*(t)$.
All of this is nice and the procedure is identical to what we have done in case of two distinct real eigenvalues a couple of lectures ago. The difference is that this general solution we constructed is complex-valued.
Question: How can we obtain a general solution that yields only (and all of the possible) real valued solutions of this system of ODEs?
Just looking at the ODE, we see that real-valued initial data will result in real valued solutions since $\mathbf{A}$ is real valued. We proceed with a key observation.
Recall that given a complex number $z= a+ i b$, where $a$ and $b$ are real numbers, $a$ is called the real part of $z$, denoted by $\text{Re}\,z$ and $b$ is called the imaginary part of $z$, denoted by $\text{Im}\,z$. Now, clearly we can express real and imaginary parts of $z$ as linear combinations of $z$ and $z^*$ as follows:
$$ \frac{1}{2}(z+z^*) = a = \text{Re}\,z $$and
$$ \frac{1}{2i}(z-z^*) = b= \text{Im}\,z $$Following this idea, we express the real and imaginary parts of the solution $\mathbf{x}_1(t)$ we obtained above in the same way. Let $\mathbf{v}(t):=\text{Re}\,\mathbf{x}_1(t)$ and $\mathbf{w}:=\text{Im}\,\mathbf{x}_1(t)$ so that
$$ \mathbf{x}_1(t)=\mathbf{v}(t) + i \mathbf{w}(t) $$Then, we mimic the idea above to get
$$ \mathbf{v}(t) = \frac{1}{2}(\mathbf{x}_1(t)+\mathbf{x}_1^*(t) ) = \frac{1}{2}\mathbf{x}_1(t) + \frac{1}{2}\mathbf{x}_2(t), $$and
$$ \mathbf{w}(t) = \frac{1}{2i}(\mathbf{x}_1(t)-\mathbf{x}_1^*(t) ) = \frac{1}{2i}\mathbf{x}_1(t) - \frac{1}{2i}\mathbf{x}_2(t). $$Now, observe that $\mathbf{v}(t)$ is simply obtained from the complex general solution by setting the constants $c_1=c_2=\frac{1}{2}$ and $\mathbf{w}(t)$ is simply obtained from the complex general solution by setting the constants $c_1=\frac{1}{2i}$, $c_2=\frac{-1}{2i}$. Bot both of these vector solutions are real-valued and their Wronskian is nonzero. We arrived at the follwoing fact:
Fact: The real and imaginary parts of $\mathbf{x}_1(t)$ ($\mathbf{v}(t)$ and $\mathbf{w}(t)$, respectively) are independent real-valued solutions of the system of ODEs. Therefore they can be used as building blocks of the real-valued general solution:
$$ \mathbf{x}(t) = c_1\mathbf{v}(t) + c_2\mathbf{w}(t). $$Now, to find the general solution of the ODE at the beginning, we need to calculate the real and imaginary parts of
$$ \mathbf{x}_{1}(t) = e^{\lambda_1 t}\mathbf{u}_1 = e^{(1+2i)t}\begin{bmatrix}1 \\ 1-i\end{bmatrix}. $$Recall the De Moivre's formula: $e^{i\theta}=\cos\theta + i \sin\theta$.
Then
$$ \mathbf{x}_{1}(t)=e^{(1+2i)t}\begin{bmatrix}1 \\ 1-i\end{bmatrix}=e^t(\cos(2t)+i\sin(2t))\begin{bmatrix}1 \\ 1-i\end{bmatrix}=e^t\begin{bmatrix}\cos(2t)+i\sin(2t) \\ \cos(2t)+i\sin(2t)-i\cos(2t)+\sin(2t)\end{bmatrix} $$which can be written as
$$ \mathbf{x}_1(t)=\underbrace{e^t\begin{bmatrix}\cos(2t) \\ \cos(2t) + \sin(2t)\end{bmatrix}}_{\text{$\mathbf{v}(t)$}}+i\cdot \underbrace{e^t\begin{bmatrix}\sin(2t) \\ \sin(2t) - \cos(2t)\end{bmatrix}}_{\text{$\mathbf{w}(t)$}} $$Taking the linear combination of the real and imaginary parts above, we obtain the real-valued general solution:
$$ \mathbf{x}(t)=c_1e^t\begin{bmatrix}\cos(2t) \\ \cos(2t) + \sin(2t)\end{bmatrix}+c_2 e^t\begin{bmatrix}\sin(2t) \\ \sin(2t) - \cos(2t)\end{bmatrix}. $$Note that all solutions (for all choices of $c_1$ and $c_2$) grow exponentially as $t\to+\infty$. This behavior is solely governed by the sign of the real part of the eigenvalues found. Since $\lambda_1 =1+2i$ and $\lambda_2 =1-2i$, $\text{Re}\,\lambda_1=\text{Re}\,\lambda_2=1$, which is the coefficient of $t$ in the exponential term above. If the real parts of the eigenvalues were negative, then we would have solutions that decay exponentially in time. In any case, the solutions exhibit oscillatory behavior due to the presence of sine and cosine functions.
You can find below the phase plane. Due to the presence of oscillatory sine and cosine terms, the trajectories are spirals. The origin (which is the only equilibrium solution) is called the spiral source in this case. The direction of the spiral (clockwise vs. counter-clockwise) can be determined by picking a particular solution (choosing some values for $c_1$ and $c_2$) and evaluating it at $t=0$. Once $\mathbf{x}(0)$ is obtained, one can find the direction vector at that point in the phase plane: $\mathbf{x}'(0)=\mathbf{A}\mathbf{x}(0)$. Please see the hand-written notes 'Here is how I would do it' that I sent.
In [6]:
using PlotlyJS
tt1 = -5:0.01:3;
tt = 0:0.1:5;
v1 = exp.(tt1).*cos.(2*tt1);
v2 = exp.(tt1).*(cos.(2*tt1)+sin.(2*tt1));
w1 = exp.(tt1).*sin.(2*tt1);
w2 = exp.(tt1).*(sin.(2*tt1)-cos.(2*tt1));
tracesODE1 = GenericTrace[];
trace = scatter(x=v1, y=v2, name = "v(t)");
push!(tracesODE1, trace);
# trace = scatter(x=-v1, y=-v2, name = "v(t)");
# push!(tracesODE1, trace);
trace = scatter(x=w1, y=w2, name = "w(t)");
push!(tracesODE1, trace);
# trace = scatter(x=-w1, y=-w2, name = "w(t)");
# push!(tracesODE1, trace);
c1=1;
c2=-1/2;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=x13, y=x23, name = "a solution");
push!(tracesODE1, trace);
c1=0.6;
c2=2;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=x13, y=x23, name = "another solution");
push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Trajectories in the phase plane. All solutions grow exponentially.", xaxis=attr(title="x<sub>1</sub>"),yaxis=attr(title="x<sub>2</sub>")))
Out[6]:
Here are the component plots of a solution with $c_1=-1$ and $c_2=5$.
In [5]:
using PlotlyJS
tt = 0:0.01:6;
v1 = exp.(tt1).*cos.(2*tt1);
v2 = exp.(tt1).*(cos.(2*tt1)+sin.(2*tt1));
w1 = exp.(tt1).*sin.(2*tt1);
w2 = exp.(tt1).*(sin.(2*tt1)-cos.(2*tt1));
tracesODE1 = GenericTrace[];
c1=-1;
c2=5;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=tt, y=x13, name = "x<sub>1</sub>(t)");
push!(tracesODE1, trace);
trace = scatter(x=tt, y=x23, name = "x<sub>2</sub>(t)");
push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Component plot of an oscillatory and exponentially growing solution.", xaxis=attr(title="t"),yaxis=attr(title="x<sub>1</sub>,x<sub>2</sub>")))
Out[5]:
We now consider the case where the eigenvalues are purely imaginary: their real parts are zero. In this case, there is no exponential time-dependence and we saw that the solution trajectories are confined to elliptical orbits. Here is such a system:
$$ \mathbf{x}'=\begin{bmatrix}2 & -5 \\ 1 & -2\end{bmatrix}\mathbf{x}. $$The eigenvalues of the coefficient matrix are $\lambda_1=i$ and $\lambda_2=-i$. The eigenvector for $\lambda_1=i$ is
$$ \mathbf{u}_1=\begin{bmatrix}2+i \\ 1\end{bmatrix}, $$and then a solution of the system of ODEs is
$$ \mathbf{x}_1(t)=e^{i t}\begin{bmatrix}2+i \\ 1\end{bmatrix}=\underbrace{\begin{bmatrix}2\cos t - \sin t\\ \cos t \end{bmatrix}}_{\text{$\mathbf{v}(t)$}}+ i\underbrace{\begin{bmatrix}2\sin t + \cos t\\ \sin t \end{bmatrix}}_{\text{$\mathbf{w}(t)$}}. $$As before, from the real and imaginary parts of this solution we can build a real-valued general solution:
$$ \mathbf{x}(t)=c_1\mathbf{v}(t)+c_2 \mathbf{w}(t)=c_1\begin{bmatrix}2\cos t - \sin t\\ \cos t \end{bmatrix}+c_2\begin{bmatrix}2\sin t + \cos t\\ \sin t \end{bmatrix} $$Observe that for any choice of $c_1$ and $c_2$ both components of the solution are periodic functions of $t$. Therefore the trajectory of any solution is an ellipse in the phase plane. Below is a phase plane portrait. The equilibrium is called a stable center.
In [8]:
using PlotlyJS
tt = 0:0.01:10;
v1 = 2.*cos.(tt) - sin.(tt);
v2 = cos.(tt);
w1 = 2.*sin.(tt) + cos.(tt);
w2 = sin.(tt);;
tracesODE1 = GenericTrace[];
trace = scatter(x=v1, y=v2, name = "v(t)");
push!(tracesODE1, trace);
trace = scatter(x=w1, y=w2, name = "w(t)");
push!(tracesODE1, trace);
c1=1/3;
c2=-1/2;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=x13, y=x23, name = "a solution");
push!(tracesODE1, trace);
c1=0.6;
c2=2;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=x13, y=x23, name = "another solution");
push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Trajectories in the phase plane. All solutions have elliptical orbits.", xaxis=attr(title="x<sub>1</sub>"),yaxis=attr(title="x<sub>2</sub>")))
Out[8]:
The component plot of a solution with $c_1=2$ and $c_2=-1/3$ is below.
In [9]:
using PlotlyJS
tt = 0:0.01:20;
v1 = 2.*cos.(tt) - sin.(tt);
v2 = cos.(tt);
w1 = 2.*sin.(tt) + cos.(tt);
w2 = sin.(tt);;
tracesODE1 = GenericTrace[];
c1=2.0;
c2=-1/3;
x13 =c1*v1+c2*w1;
x23 =c1*v2+c2*w2;
trace = scatter(x=tt, y=x13, name = "x<sub>1</sub>(t)");
push!(tracesODE1, trace);
trace = scatter(x=tt, y=x23, name = "x<sub>2</sub>(t)");
push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Component plot of an oscillatory and exponentially growing solution.", xaxis=attr(title="t"),yaxis=attr(title="x<sub>1</sub>,x<sub>2</sub>")))
Out[9]: