There are often scenarios where we are interested in understanding how two quantities $y_1$ and $y_2$ that change simultaneously evolve in time. For example, if the system of equations that model this change is
\begin{align} \frac{dy_1}{dt} &= -2 y_1\\ \frac{dy_2}{dt} &= 5 y_2, \end{align}the equations are decoupled: $y_1$ does not communicate with $y_2$. We can integrate immediately and get
\begin{align} y_1(t) &= C_1 e^{-2 t}\\ y_2(t) &= C_2 e^{5 t}. \end{align}But what if the rate of change of $y_1$ dependent on the current value of $y_2$ instead of depending only on the value of $y_1$? For example:
\begin{align} \frac{dy_1}{dt} &= -2 y_1 + 3 y_2 - 12\\ \frac{dy_2}{dt} &= y_1 + 5 y_2. \end{align}We will cover how can we solve this system of first order linear differential equations in the next couple of lectures. For now we introduce some terminology.
Observe that we can write this system as a matrix equation:
$$ \begin{bmatrix}\frac{dy_1}{dt}\\\frac{dy_2}{dt}\end{bmatrix}= \begin{bmatrix} -2 & 3 \\ 1 & 5 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} + \begin{bmatrix} -12 \\ 0 \end{bmatrix}. $$The matrix that consists of the coefficients of the unknowns $y_1$ and $y_2$ on the right hand side of the system of ODEs is called the coefficient matrix. In this case this matrix is constant, it does not depend on the independent variable $t$. But in general, it might as well depend on $t$, such as:
\begin{align} \begin{bmatrix}\frac{dy_1}{dt}\\ \frac{dy_2}{dt}\end{bmatrix}&= \begin{bmatrix} p_{11}(t)y_1(t) + p_{12}(t)y_2(t) + g_1(t) \\ p_{21}(t)y_1(t) + p_{22}(t)y_2(t) + g_2(t) \end{bmatrix}\\ &=\mathbf{P}(t)\mathbf{y}(t) + \mathbf{g}(t),~\mathbf{y}(t_0)=\mathbf{y}_0. \end{align}Theorem: If all of the functions $p_{ij}(t)$ and $g_{j}(t)$ are continuous in an open interval $I$ that contains the initial point $t_0$, then the IVP above has a unique solution on the interval $I$.
Given a system of first order differential equations like above:
$$ \mathbf{u}'(t) = \mathbf{K}\mathbf{u}(t) + \mathbf{b}, $$the components $u_1$, $u_2$ of the solution $\mathbf{u}$ are called state variables since the solution $\mathbf{u}$ describes the state of the system.
The $(u_1,u_2)$-plane is called the phase plane (or phase space). Given an initial condition
$$ \mathbf{u}_0 = \begin{bmatrix}u_1(0) \\ u_2(0) \end{bmatrix}, $$the corresponding solution leaves a path
$$ \begin{bmatrix}u_1(t) \\ u_2(t) \end{bmatrix} $$on the phase plane. Such paths are called solution trajectories or orbits.
$$ \begin{bmatrix}u_1(t) \\ u_2(t) \end{bmatrix}=u_1(t)\mathbf{i} + u_2(t)\mathbf{j} $$can also be interpreted as a position vector of a particle moving in the phase space according to the dynamics of the ODE system.
The plots of the components $u_1(t)$ and $u_2(t)$ of the solution vs. $t$ are called component plots.
Consider the initial value problem
\begin{align} u_1' &= u_1 - 4u_2 + 2t,~u_1(0)=1\\ u_2' &= u_1 - 3u_2 + -3,~u_2(0)=-2. \end{align}It can be written as
$$ \mathbf{u}'(t)=\begin{bmatrix} 1&-4\\1&-3\end{bmatrix}\mathbf{u}+\begin{bmatrix}2t\\-3\end{bmatrix},~\mathbf{u}(0)=\begin{bmatrix}1\\-2 \end{bmatrix}, $$the coefficient matrix is
$$ \mathbf{K} = \begin{bmatrix} 1&-4\\1&-3\end{bmatrix}. $$Given that the solution of the IVP is
$$ \mathbf{u}(t)=e^{-t}\begin{bmatrix}2t-1\\t-1\end{bmatrix}+\begin{bmatrix}6t+2\\2t-1\end{bmatrix}, $$we can produce a component plot:
In [1]:
using PlotlyJS;
# using Plots
# plotly()
tt = 0:0.001:12;
tracesODE1 = GenericTrace[];
u10=1.0
u20=-2.0
trace = scatter(x=tt, y=exp(-tt).*(2.0*tt-1.0)+6.0*tt+2.0, name = "u<sub>1</sub>(0)=$(u10)");
push!(tracesODE1, trace);
trace = scatter(x=tt, y=exp(-tt).*(tt-1.0)+2.0*tt-1.0, name = "u<sub>2</sub>(0)=$(u20)");
push!(tracesODE1, trace);
trace=scatter(x=[0],y=[1.0], marker_size=12, name="u<sub>1</sub>(0)");
push!(tracesODE1, trace);
trace=scatter(x=[0],y=[-2.0], marker_size=12, name="u<sub>2</sub>(0)");
push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Component plots", xaxis=attr(title="t"),yaxis=attr(title="u1,u2")))
Out[1]:
And the solution trajectory:
In [2]:
tt = 0:0.001:12;
u1 = exp(-tt).*(2.0*tt-1.0)+6.0*tt+2.0;
u2 = exp(-tt).*(tt-1.0)+2.0*tt-1.0;
tracesODE1 = GenericTrace[];
u10=1.0
u20=-2.0
trace = scatter(x=u1, y=u2, name = "(u1,u2)");
push!(tracesODE1, trace);
# trace = scatter(x=tt, y=exp(-tt).*(tt-1.0)+2.0*tt-1.0, name = "u<sub>2</sub>(0)=$(u20)");
# push!(tracesODE1, trace);
trace=scatter(x=[1],y=[-2.0], marker_size=12, name="Initial value");
push!(tracesODE1, trace);
# trace=scatter(x=[0],y=[-2.0], marker_size=12, name="u<sub>2</sub>(0)");
# push!(tracesODE1, trace);
plot(tracesODE1, Layout(title="Orbit (trajectory) in the phase plane", xaxis=attr(title="u<sub>1</sub>"),yaxis=attr(title="u<sub>2</sub>")))
Out[2]:
Now we consider systems of the form
$$ \mathbf{x}'(t) = \mathbf{A}\,\mathbf{x}(t)+\mathbf{b}. $$Note that by the existence and uniqueness theorem we covered, the corresponding initial value problem has a unique solution $-\infty \lt t \lt +\infty$. Just as we did in case of a single autonomous differential equation, we can look for equilibrium solutions, i.e. solutions that do not change in time. Thus we seek solutions of the above system that satisfy
$$ \mathbf{x}'(t) =\begin{bmatrix}0\\0\end{bmatrix}, $$which is possible if and only if
$$ \mathbf{A}\,\mathbf{x}(t) = -\mathbf{b}. $$Note that this equation has a unique solution if and only if $\det(\mathbf{A})\neq 0$, i.e. the coefficient matrix is nonsingular. We arrived at the following fact.
Fact: $\mathbf{x}'(t) = \mathbf{A}\,\mathbf{x}(t)+\mathbf{b}$ has a unique equilibrium solution $\mathbf{x}_{\text{eq}}$ if and only if $\det(\mathbf{A})\neq 0$. In this case
$$ \mathbf{x}_{\text{eq}}=-\mathbf{A}^{-1}\mathbf{b} $$On the other hand, if $\det(\mathbf{A})=0$, either there is no equilibrium solution or there are infinitely many of them.
It is important to observe that to find equilibrium solutions of a linear autonomous system, we solve an algebraic equation. We will come back to solving linear autonomous systems later.
Consider the IVP
$$ \mathbf{u}' = \begin{bmatrix}\frac{-13}{8}&\frac{3}{4}\\\frac{1}{4}&\frac{-1}{4}\end{bmatrix}\mathbf{u} + \begin{bmatrix} 14 \\ 0 \end{bmatrix}, $$$$ \mathbf{u}(0) = \begin{bmatrix} 36 \\ 0 \end{bmatrix}. $$Note that the coefficient matrix
$$ \mathbf{A}=\begin{bmatrix}\frac{-13}{8}&\frac{3}{4}\\\frac{1}{4}&\frac{-1}{4}\end{bmatrix} $$is nonsingular since
$$ \det(\mathbf{A})=\frac{7}{32}\neq 0. $$Therefore this systems has exactly one equilibrium solution given by
$$ \mathbf{x}_{\text{eq}}=-\mathbf{A}^{-1}\begin{bmatrix} 14 \\ 0 \end{bmatrix}. $$We proceed with calculating the inverse of the coefficient matrix:
$$ \mathbf{A}^{-1} = \frac{32}{7}\begin{bmatrix}\frac{-1}{4}&\frac{-3}{4}\\\frac{-1}{4}&\frac{-13}{8}\end{bmatrix}. $$Then
$$ \mathbf{x}_{\text{eq}} = -\frac{32}{7}\begin{bmatrix}\frac{-1}{4}&\frac{-3}{4}\\\frac{-1}{4}&\frac{-13}{8}\end{bmatrix}\begin{bmatrix} 14 \\ 0 \end{bmatrix}=\begin{bmatrix} 16 \\ 16 \end{bmatrix}. $$Consider the initial value problem
\begin{align} y'' + p(t)y' + q(t)y &= g(t)\\ y(0)&=y_0\\ y'(0)&=y_1. \end{align}Set $x_1=y$ and $x_2=y'$ and write a system fo the unknown vector $\mathbf{x}=\begin{bmatrix}x_1\\x_2 \end{bmatrix}$. Note that the derivative of this unknown vector contains the term $y''$:
$$ \frac{d}{dt}\begin{bmatrix}x_1\\x_2 \end{bmatrix}=\begin{bmatrix}y'\\y'' \end{bmatrix}. $$and the initial data can be written as
$$ \begin{bmatrix}x_1(0)\\x_2(0) \end{bmatrix}=\begin{bmatrix}y_0\\y_1 \end{bmatrix}. $$Then we can write
$$ \mathbf{x}' = \begin{bmatrix}0 & 1 \\ -q(t) & -p(t) \end{bmatrix}\mathbf{x} + \begin{bmatrix}0\\g(t) \end{bmatrix},~\mathbf{x}(0)=\begin{bmatrix}y_0\\y_1 \end{bmatrix}. $$If we solve this system, we recover the solution of the IVP problem for the second order ODE from the first component of the solution of the system.
Exercise: Everybody should check that the system above is equivalent to the second order ODE.
In [ ]: