3  Systems of First order equations

3.1 The General Nonlinear System: Existence and Uniqueness of Solutions

Consider a system of differential equations of the form \(y'(t) = f(t,y(t))\), where \(f:{\mathbb R}\times{\mathbb R}^n\to{\mathbb R}^n\) and we seek a solution \(y:{\mathbb R}\to {\mathbb R}^n\). That is, \(y\) is a vector, or n-tuple of functions, and the function \(f\) takes a vector and a real number as inputs and returns a vector as an output. The equations are often accompanied by initial conditions of the form \(y(t_0) = y_0\). The solution is a curve in \({\mathbb R}^n\) parameterized by \(t\) and passing through the point \(y_0\).

Sufficient conditions to guarantee the existence of a unique solution through any point are that \(f\) be continuous in \(t\) and differentiable in \(y\). However, a weaker condition than differentiability is also known to suffice.

Definition 3.1 (Lipschitz Continuity) A function \(f\) is said to be Lipschitz continuous on some subset \(U\) of its domain if there exists a constant \(M\) for which \(\left\|f(x)-f(y)\right\| < M\left\|x-y\right\|\) whenever \(x,y\in U\) with \(x\ne y\).

Equivalently, a function is Lipschitz continuous if the slopes of its secant lines are all bounded. For example, \(|x|\) is Lipschitz continuous; \(x^{2/3}\) is continuous but not Lipschitz since the secant lines tend towards vertical near the origin. Note, \(x^{2/3}\) is (locally) Lipschitz so long as the set of interest is bounded away from the origin. In one dimension, \(\left\|x-y\right\|\) can be assumed to be \(|x-y|\). In higher dimensions \(\left\|\cdot\right\|\) is most often the standard Euclidian distance.

The classic result for the existence of solutions to our system of differential equations is as follows.

Theorem 3.1 If \(f(t,y)\) is continuous in the independent variable, \(t\), and Lipschitz continuous in \(y\), then the initial value problem \(y' = f(t,y)\), \(y(t_0) = y_0\) has a unique solution on some open interval \(a < t_0 < b\).

Various proofs of the above theorem can be found in most advanced books on differential equations.

Example 3.1 An example for which solutions exist, but are not unique is the simple scalar initial value problem \(y' = y^{2/3}\), \(y(0) = 0\). The function \(y^{2/3}\) is continuous for all \(y\), but has a cusp at the origin. It is not differentiable, nor is it Lipschitz continuous at the origin. Yet the differential equation is easily solved using the method of separation of variables, leading to the solutions \(y(t) = \frac{1}{3}(t+c)^3\). Taking \(c=0\) gives a solution passing through the origin. A second obvious solution is \(y(t) = 0\). Moreover, we can splice these solutions together to get an infinite family of solutions to the initial value problem: \[\begin{equation} y(t) = \begin{cases} \frac{1}{3}(t+c_1)^3 & t < -c_1 \\ 0 & -c_1 \le t \le c_2 \\ \frac{1}{3}(t-c_2)^3 & t > c_2 \end{cases} \end{equation}\] with \(c_1\) and \(c_2\) arbitrary nonnegative constants.

Systems of higher order differential equations can be expressed as systems of linear equations by introducing the derivatives as new state variables.

Example 3.2 Express the second order equation \((t-3)y''(t) + t^2y'(t) + 7y(t) = f(t)\) as a system of two first order equations.

Let \(z = y'\). Then the second order equation becomes \((t-3)z'(t) + t^2z(t) + 7y(t) = f(t)\). This together with the definition of \(z\) gives a pair of first order equations. Dividing through by \((t-3)\) puts the system in the basic form \[\begin{align*} y'(t) &= z(t) \\ z'(t) &= -\frac{t^2}{t-3}z(t) - \frac{7}{t-3}y(t) + f(t) \end{align*}\] The right hand side is linear, and hence differentiable, in \((y,z)\) and continuous in \(t\) on any interval bounded away from \(t=3\).

4 Homogeneous Constant Coefficient Linear Systems

4.1 The big picture

Consider the linear system of ordinary differential equations \(x' = Ax\) where \(x(t) = \left( x_1(t), x_2(t), \dots, x_n(t) \right)^T\)and \(A\) is an \(n\times n\) matrix of real numbers (constants). The general theory for finding all solutions to the linear system is based on an eigenvalue decomposition of the matrix \(A\). Let \(r_1, r_2, \dots r_m\) be the \(m\) eigenvalues of \(A\). The theory of matrices and eigenvalues guarantees that \(A\) has at least one eigenvalue, and no more than \(n\) distinct eigenvalues, hence we can arrange things so that \(1 \le m \le n\) and all the \(m\) eigenvalues are distinct. That is, \(r_i=r_j \leftrightarrow i=j\). A handy theorem from linear algebra states that there is an invertible matrix \(S\), for which \(S^{-1}AS = \begin{pmatrix} J_1 & & 0 \\ & \ddots & \\ 0 & & J_m \end{pmatrix}\).

We know from linear algebra that each eigenvalue of \(A\) has a multiplicity.

4.2 Real and Distinct Eigenvalues

Consider the system \(y'=Ay\) where \(A\) is a real \(n\times n\) matrix, and suppose that \(A\) has a real eigenvalue \(r\) with an associated eigenvector \(u\). Then \(y(t) = e^{rt}u\) a solution. This is easy to verify by direct substitution.

Suppose \(u_1,\dots,u_m\) are \(m\) linearly independent eigenvectors of \(A\) with associated eigenvalues \(r_1,\dots,r_m\). We look for a solution of the form \[\begin{equation}y(t) = \sum_i^m x_i(t)u_i\end{equation}\] Substituting this form into the differential equation leads to \[\begin{equation} \sum_i^m x'_i(t)u_i = A\sum_i^m x_i(t)u_i = \sum_i^m r_ix_i(t)u_i \Rightarrow \sum_i^m \left(x'_i(t)-rx_i(t)\right)u_i = 0 \end{equation}\] Since the vectors are linearly independent, it follows that \(x'_i(t)-rx_i(t) = 0\), \(i=1,\dots,m\), which has the solutions \(x_i(t) = c_ie^{r_it}\), \(i=1,\dots,m\), where each \(c_i\) is an arbitrary constant. Thus, \[\begin{equation}y(t) = \sum_i^m c_iu_ie^{r_it} \end{equation}\] The effect of choosing the eigenvectors as a basis is to transform the system to a set of \(m\) decoupled equations that have well-known solutions.

4.3 Complex Eigenvalues

Consider the system \(y'(t) = Ay(t)\) where \(A\) is a real \(n\times n\) matrix, and suppose that \(A\) has a complex eigenvalue \(r = a+bi\) with an associated eigenvector \(u = v + iw\). It follows by direct substitution that \(v-iw\) is also an eigenvector of \(A\) associated with the eigenvalue \(a-ib\).

We seek a solution of the form \(y(t) = x_1(t) v + x_2(t) w\) where \(x_1\) and \(x_2\) are scalar functions of \(t\). Direct computation yields that \(A(x_1v + x_2w) = (ax_1 +bx_2)v + (-bx_1+ax_2)w\). Hence, \(x_1\) and \(x_2\) must satisfy the pair of differential equations \[\begin{align*} x_1'(t) &= ax_1(t) +bx_2(t) \\ x_2'(t) &= -bx_1(t) +ax_2(t) \end{align*}\] Since \(e^{at}\) is an integrating factor for both the above equations, \[\begin{align} \frac{d}{dt} \left(x_1e^{-at}\right) &= bx_2e^{-at}, \label{intfactx1}\\ \frac{d}{dt} \left(x_2e^{-at}\right) &= -bx_1e^{-at}.\label{intfactx2} \end{align}\] Differentiating the first of these equations and eliminating \(x_2\) using the second yields \[\begin{equation} \frac{d^2}{dt^2} \left(x_1e^{-at}\right) = b\frac{d}{dt}\left( x_2e^{-at}\right) = -b^2x_1e^{-at}. \end{equation}\] Hence, \(x_1e^{-at} = c_1\cos(bt) + c_2\sin(bt)\), which implies \(x_1(t)= e^{at}\left( c_1\cos(bt) + c_2\sin(bt)\right)\). To compute \(x_2\) we use as follows: \[\begin{align*} bx_2e^{-at} &= \frac{d}{dt} \left(x_1e^{-at}\right) \\ &= \frac{d}{dt} \left( c_1\cos(bt) + c_2\sin(bt)\right) \\ &= \left( -bc_1\sin(bt) + bc_2\cos(bt)\right) \\ x_2(t) &= e^{at} \left( -c_1\sin(bt) + c_2\cos(bt)\right) \end{align*}\] A nice way to memorize these solutions is to write them in a matrix form. \[\begin{equation} \begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} = e^{at}\begin{pmatrix} \cos(bt) & \sin(bt) \\ -\sin(bt) & \cos(bt) \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} \end{equation}\] The matrix on the right is the rotation matrix corresponding to counterclockwise rotation by the angle \(bt\). Hence, the solutions form a spiral in the \(u\)-\(v\) plane.

Returning to the original coordinates, \[\begin{equation} y(t) = c_1e^{at}\left(v\cos(bt) - w\sin(bt)\right) + c_2e^{at}\left(v\sin(bt) + w\cos(bt)\right) \end{equation}\] Note that up until now we have not specified the size of the original system. The solutions found above work regardless of the number of equations in the original system. If \(A\) is a \(2\times2\) matrix, then the solution can be expressed compactly by defining \(S\) as the matrix whose columns are \(v\) and \(w\). In this case \(y = Sx = e^{at}SR(bt)C\) where \(R(\theta)\) is the rotation matrix through angle \(\theta\) and \(C\) is the transpose of \(\begin{pmatrix} c_1 & c_2\end{pmatrix}\). If we have initial conditions at \(t=0\), then \(y(0) = SC\), since \(R(0)\) is the identity (rotation by zero). Thus, \(C= S^{-1}y(0)\) and \[\begin{equation} y(t) = e^{at}SR(bt)S^{-1}y(0). \end{equation}\]

4.4 Repeated Eigenvalues

We have seen how to find solutions if we know a sufficient number of linearly independent eigenvectors. In this section we tackle the case where the geometric multiplicity of an eigenvalue (the dimension of the Eigenspace) is less than its algebraic multiplicity. We do this by first finding a general solution for a simple case and then derive a method to transform, through a change of basis, any other system to a similar simple system.

4.4.1 The second simplest case

Consider the system \[\begin{align} \dot{x}_1(t) &= rx_1(t) + x_2(t), \label{simple1}\\ \dot{x}_2(t) &= rx_2(t), \label{simple2} \end{align}\] which when written in matrix form looks like \[ \begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \end{pmatrix} = \begin{pmatrix} r & 1 \\ 0 & r \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \]

The matrix \(J = \begin{pmatrix} r & 1 \\ 0 & r \end{pmatrix}\) has a single eigenvalue \(r\) with multiplicity 2. However, all eigenvectors of \(J\) are multiples of \(\vec{u} = \begin{pmatrix}1 \\\ 0 \end{pmatrix}\). Hence, the eigenvalue \(r\) has geometric multiplicity one.

To solve the system, first integrate to obtain \(x_2(t) = c_2e^{rt}\), then substitute this into and integrate a second time using the integrating factor \(e^{-rt}\) to obtain \[\begin{align*} \dot{x}_1(t) = rx_1(t) &+ c_2e^{rt} \\ \dot{x}_1(t) - rx_1(t) &= c_2e^{rt} \\ \frac{d}{dt}\left({x}_1(t)e^{-rt}\right) &= c_2 \\ {x}_1(t)e^{-rt} &= c_2t + c_1\\ {x}_1(t) &= c_2te^{rt} + c_1e^{rt} \end{align*}\] This solution can be written in vector form as \[\begin{equation} \begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} = c_1 e^{rt}\begin{pmatrix} 1 \\ 0 \end{pmatrix} + c_2 \left( te^{rt}\begin{pmatrix} 1 \\ 0 \end{pmatrix} + e^{rt}\begin{pmatrix} 0 \\ 1 \end{pmatrix}\right) \end{equation}\] or in matrix form as \[\begin{equation} \begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} = \begin{pmatrix} e^{rt} & te^{rt}\\ 0 & e^{rt} \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} \end{equation}\] Both forms will be convenient at times.

4.4.2 Generalized Eigenvectors

Definition 4.1 (Generalized eigenvector) A generalized eigenvector of \(A\) is a vector \(v\) for which \[\left(A-rI\right)^m v = 0\] for some positive integer \(m\).

Comments:

• Clearly, any eigenvector of \(A\) is a generalized eigenvector (with \(m=1\)).
• If \(v\) is a generalized eigenvector of \(A\), then the scalar \(\lambda\) in the definition is an eigenvalue of \(A\). It also follows that \(u = \left(A-\lambda I\right)^{m-1}v\) is an associated eigenvector.

Theorem 4.1 If \(\lambda\) is an eigenvalue of \(A\) with algebraic multiplicity \(k\), then the null space of \((A-\lambda I)^k\) has dimension \(k\).

Theorem 4.2 Every matrix \(A\) has a full set of generalized eigenvectors. That is, there is a (possibly complex) linearly independent set \(v_1,\dots,v_n\) of generalized eigenvectors for every \(n\times n\) matrix \(A\).

Example 4.1 Consider the linear system \(\vec{y}' = A\vec{y}\) with \(\vec{y}(0) = \vec{y}_0\) and \[\begin{equation} A = \begin{pmatrix} 4 & -4 \\ 1 & 0 \end{pmatrix} \end{equation}\] The characteristic polynomial for \(A\) is \((4-r)(-r) -(-4)(1) = r^2 - 4r + 4 = (r-2)^2\). Hence, \(r_1 = 2\) is an eigenvalue of \(A\) with algebraic multiplicity 2. To find the eigenvectors of \(A\) associated with \(r_1\) we construct and simplify

\[(A -r_1I) = \begin{pmatrix} 4-2 & -4 \\ 1 & 0-2 \end{pmatrix} = \begin{pmatrix} 2 & -4 \\ 1 & -2 \end{pmatrix} \sim \begin{pmatrix} 1 & -2 \\ 0 & 0 \end{pmatrix}\] Thus the solutions to \((A-2I)\vec{u} = 0\) are all multiples of \(\vec{u} = \begin{pmatrix} 2 & 1 \end{pmatrix}.\) This is a one dimensional subspace. Hence \(r_1\) has geometric multiplicity 1. As it turns out, the particular generalized eigenvector we want is a solution to \[(A -r_1I)\vec{v} = \vec{u}.\] \(\vec{v}\) will be a generalized eigenvector since \[(A -r_1I)^2\vec{v} = (A -r_1I)\vec{u} = 0.\] To find \(\vec{v}\), we reduce the augmented matrix \[\left(\begin{array}{cc|c} 2 & -4 & 2 \\ 1 & -2 & 1 \end{array} \right) \sim \left(\begin{array}{cc|c} 1 & -2 & 1\\ 0 & 0 & 0 \end{array} \right)\] Thus, the entries of \(\vec{v}\) satisfy \(v_1-2v_2 = 1\). Setting \(v_2 = a \in {\mathbb R}\) give \(\vec{v} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \begin{pmatrix} 2 \\ 1 \end{pmatrix}a\). We will choose \(a=0\) to give \(\vec{v} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}\). Note that \(\vec{u}\) and \(\vec{v}\) are linearly independent and hence form a basis for \({\mathbb R}^2\).

We look for solutions relative to this basis: \[\vec{y}(t) = x_1(t)\vec{u} + x_2(t)\vec{v}\] That is, \(x_1\) and \(x_2\) are the coordinates of our solution with respect to the basis \(\left\{ \vec{u},\vec{v}\right\}\). Substituting this form for \(\vec{y}\) into the differential equation leads to \[\begin{align*} x'_1(t)\vec{u} + x'_2(t)\vec{v} &= x_1(t)A\vec{u} + x_2(t)A\vec{v} \\ &= x_1(t)r_1\vec{u} + x_2(t)\left(r_1\vec{v} + \vec{u}\right) \\ &= \left(r_1x_1(t) + x_2(t)\right)\vec{u} + r_1x_2(t)\vec{v} \end{align*}\] Finally, equating the coefficients on both sides of the equation leads to the simpler system given by Equations (\(\ref{simple1}\)) and (\(\ref{simple2}\)) above, which have solutions \[\begin{align*} x_1(t) &= c_2te^{r_1t} + c_1e^{r_1t},\\ x_2(t) &= c_2e^{r_1t}. \end{align*}\] This leads to a solution for \(\vec{y}\) of \[\vec{y}(t) = \left(c_2te^{r_1t} + c_1e^{r_1t}\right)\vec{u} + c_2e^{r_1t}\vec{v} = c_1e^{r_1t}\vec{u} + c_2\left(te^{r_1t}\vec{u} + e^{r_1t}\vec{v}\right)\] For our example we have \(r_1 = 2\), \(\vec{u} = \begin{pmatrix} 2 \\ 1\end{pmatrix}\) and \(\vec{v} = \begin{pmatrix} 1 \\ 0\end{pmatrix}\). Hence \[\vec{y}(t) = c_1e^{r_1t}\begin{pmatrix} 2 \\ 1\end{pmatrix} + c_2\left(te^{r_1t}\begin{pmatrix} 2 \\ 1\end{pmatrix} + e^{r_1t}\begin{pmatrix} 1 \\ 0\end{pmatrix}\right)\]

4.4.3 The Jordan form

The question remains if we have found all the solutions. The answer to this lies in a theorem from linear algebra that states, roughly, that given any \(n\times n\) matrix \(A\), there is a linearly independent set of generalized eigenvectors of \(A\) that spans \({\mathbb R}^n\). Moreover, proceeding as we have above, we will find such a set. It can then be shown that the solutions we’ve found will be linearly independent and yield the complete general solution to \(y' = Ay\).

5 The nonhomogeneous problem

5.1 The method of variation of parameters

Consider the system \[\begin{equation} \frac{d}{dt}y(t) = Ay(t) + f(t), \qquad y(t_0) = y_0, \label{nonhom} \end{equation}\] where \(A\) is a given \(n\times n\) matrix of real numbers and \(f\) is a given vector-valued function from \({\mathbb R}\) to \({\mathbb R}^n\) Suppose we know \(n\) linearly independent solutions to the homogeneous system \(y' = Ay\). Then the fundamental matrix, \(\Psi\), whose columns are these \(n\) solutions satisfies \[\begin{equation} \frac{d}{dt}\Psi(t) = A\Psi(t) + f(t). \end{equation}\] Applying the method of variation of parameters, we seek a solution to of the form \(y(t) = \Psi(t)u(t)\). Substituting this form into leads to \[\Psi'u + \Psi u' = A\Psi u + f.\] Since \(\Psi\) is a solution to the homogeneous problem, the above equation reduces to \[ \Psi u' = f,\] and hence, \[\begin{equation}u(t) = \int \Psi^{-1}(t) f(t) \, dt.\end{equation}\] To satisfy, the initial conditions, we choose the particular solution \(y_p = \Psi u\) satisfying \(y_p(t_0) = y_0\). Then \[\begin{equation}y(t) = \Psi(t)\Psi^{-1}(t_)y_0 + \Psi(t)\int_{t_0}^{t} \Psi^{-1}(\tau) f(\tau) \, d\tau.\end{equation}\]

Example 5.1 Consider the system \(y' = Ay + f\) with \(A = \begin{pmatrix} 3 & 3 \\ 1 & 5 \end{pmatrix}\) and \(f = \begin{pmatrix} t \\ e^{t} \end{pmatrix}.\) The characteristic polynomial of \(A\) is \((3-r)(5-r)-3 = (6-r)(6+r)\). The eigenvectors are \(r_1 = 6\) and \(r_2 = 2\) with associated eigenvectors \(u_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\) and \(u_2 = \begin{pmatrix} 3 \\ -1 \end{pmatrix}\). Thus two linearly independent solutions to the homogeneous problem are \(y_1(t) = e^{6t}\begin{pmatrix} 1 \\ 1 \end{pmatrix}\) and \(y_2(t) = e^{2t}\begin{pmatrix} 3 \\ -1 \end{pmatrix}\). One possible fundamental matrix is \[\begin{equation} \Psi(t) = \begin{pmatrix} e^{6t} & 3e^{2t} \\ e^{6t} & -e^{2t} \end{pmatrix} = \begin{pmatrix} 1 & 3 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} e^{6t} & 0 \\ 0 & e^{2t} \end{pmatrix} \end{equation}\] Then \(\Psi^{-1}(t) = \frac{1}{4}\begin{pmatrix} e^{-6t} & 3e^{-6t} \\ e^{-2t} & -e^{-2t} \end{pmatrix}\), and \[\begin{equation} \Psi(t)\Psi^{-1}(\tau) = \frac{1}{4}\begin{pmatrix} e^{6(t-\tau)}+3e^{2(t-\tau)} & 3e^{6(t-\tau)}-3e^{2(t-\tau)} \\ e^{6(t-\tau)}-e^{2(t-\tau)} & 3e^{6(t-\tau)}+e^{2(t-\tau)} \end{pmatrix}. \end{equation}\] The complementary solution satisfying the initial conditions is \[\begin{equation} y_c(t) = \Psi(t)\Psi^{-1}(0)y_0 = \frac{1}{4}\begin{pmatrix} e^{6t}+3e^{2t} & 3e^{6t}-3e^{2t} \\ e^{6t}-e^{2t} & 3e^{6t}+e^{2t} \end{pmatrix}y_0. \end{equation}\] The particular solution satisfying \(y_p(0)=0\) is \[\begin{align*} y_p(t) &= \int_0^t\Psi(t)\Psi^{-1}(\tau)f(\tau)\, d\tau \\ &= \frac{1}{4}\int_0^t \strut \begin{pmatrix} \tau e^{6(t-\tau)}+3\tau e^{2(t-\tau)} + 3e^{6t-5\tau}-3e^{2t-\tau} \\ \tau e^{6(t-\tau)}-\tau e^{2(t-\tau)} + 3e^{6t-5\tau}+e^{2t-\tau} \end{pmatrix}\begin{pmatrix} \tau \\ e^{\tau} \end{pmatrix} \, d\tau.\\ &= {\small \begin{pmatrix}1\\1\end{pmatrix}}\int_0^t\tfrac{\tau}{4} e^{6(t-\tau)} \,d\tau + {\small \begin{pmatrix}3\\-1\end{pmatrix}}\int_0^t \tfrac{\tau}{4} e^{2(t-\tau)} \,d\tau + {\small \begin{pmatrix}3\\3\end{pmatrix}}\int_0^t \tfrac{1}{4}e^{6t-5\tau} \,d\tau - {\small \begin{pmatrix}3\\-1\end{pmatrix}}\int_0^t \tfrac{1}{4}e^{2t-\tau} \,d\tau \\ &= {\tiny \begin{pmatrix}1\\1\end{pmatrix}} \left( \tfrac{1}{144}-\tfrac{t}{24} -\tfrac{1}{144}e^{6t} \right) + {\tiny \begin{pmatrix}3\\-1\end{pmatrix}} \left(\tfrac{1}{16} -\tfrac{t}{8} -\tfrac{1}{16}e^{2t} \right) + {\tiny \begin{pmatrix}3\\3\end{pmatrix}} \left(\tfrac{1}{20}e^{6t} - \tfrac{1}{20}e^t \right) - {\tiny \begin{pmatrix}3\\-1\end{pmatrix}} \left(\tfrac{1}{4}e^{2t} -\tfrac{1}{4}e^t \right) \\ &= \tfrac{113}{720} e^{6t}{\small \begin{pmatrix} 1 \\ 1 \end{pmatrix}} -\tfrac{3}{16}e^{2t} {\small \begin{pmatrix} 3 \\ -1 \end{pmatrix}} -\tfrac{t}{12}{\small \begin{pmatrix} 5 \\ -1 \end{pmatrix}} -\tfrac{1}{36}{\small \begin{pmatrix} 7 \\ 2 \end{pmatrix}} +\tfrac{1}{5}e^t {\small \begin{pmatrix} 3 \\ -2 \end{pmatrix}} \end{align*}\]

Note the first two terms are in fact homogeneous solutions that arise here to ensure the particular solution is zero when \(t=0\). The last three terms are of the form \(At + B + Ce^t\) which we would expect from the method of undetermined coefficients, with the caveat that \(A\), \(B\) and \(C\) are vectors.

The general solution is simply \(y(t) = y_c(t) + y_p(t)\) with the as-yet unspecified initial conditions taken as a pair of arbitrary constants.

5.2 The solution via Laplace Transforms

Applying the Laplace Transform to leads to \[\begin{equation} sY(s) - y_0 = AY(s) + F(s).\end{equation}\] Note that \(Y\) is a vector whose entries are the transforms of the corresponding entry of \(y\). Since the transform is a linear operation (an integration), the matrix is essential pulled out of the integral. Solving for \(Y\) we find \[\begin{equation}Y(s) = (sI-A)^{-1}\left(y_0 + F\right).\end{equation}\] Let \(G(s) = (sI-A)^{-1}\) and \(g = {\mathcal L}^{-1}\left\{G\right\}\). Note that \(G\) and \(g\) are both \(n\times n\) matrices. The solution to the differential equation can be represented as a convolution: \[\begin{equation}y(t) = g(t)y_0 + \int_0^t g(t-\tau)f(\tau)\, d\tau.\end{equation}\]

Example 5.2 Consider the system \(y' = Ay + f\) with \(A = \begin{pmatrix} 3 & 3 \\ 1 & 5 \end{pmatrix}\) and \(f = \begin{pmatrix} t \\ e^{t} \end{pmatrix}.\) \[\begin{align*} G(s) &= (sI-A)^{-1} \\ &= \begin{pmatrix} s-3 & -3 \\ -1 & s-5 \end{pmatrix}\\ &= \begin{pmatrix} \frac{s-5}{(s-6)(s-2)} &\ & \frac{3}{(s-6)(s-2)} \\[6pt] \frac{1}{(s-6)(s-2)} &\ & \frac{s-3}{(s-6)(s-2)} \end{pmatrix} \\[12pt] &= \frac{1}{4}\begin{pmatrix} \frac{1}{s-6} + \frac{3}{s-2} &\ & \frac{3}{s-6} - \frac{3}{s-2} \\[6pt] \frac{1}{s-6} - \frac{1}{s-2} &\ & \frac{3}{s-6} + \frac{1}{s-2} \end{pmatrix}\\[12pt] g(t) &= \frac{1}{4}\begin{pmatrix} e^{6t} + 3e^{2t} &\ & 3e^{6t} - 3e^{2t} \\[6pt] e^{6t} - e^{2t} &\ & 3e^{6t} + e^{2t} \end{pmatrix} \end{align*}\] Note that the system for this example is the same as that of the previous example, and that \(g(t-\tau) = \Psi(t)\Psi^{-1}(\tau)\) where \(\Psi\) is the fundamental matrix from the previous example. Of course, \(g(t) = \Psi(t)\Psi^{-1}(0)\), which means the rest of the computations follow the previous example exactly. The two methods are identical!

Example 5.3  

Solve the system \(x' = Ax + g\) where \(A = \begin{pmatrix} 4 & 3 \\ -1 & 0 \end{pmatrix}\) and \(g = \begin{pmatrix} 0 \\ t \end{pmatrix}\).

The matrix \(A\) has two eigenvalues, \(r_1 = 1\) and \(r_2=3\). Two associated eigenvectors are \(u_1 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}\) and \(u_2 = \begin{pmatrix} 3 \\ -1 \end{pmatrix}\), respectively. Thus, a fundamental matrix for the system is \[\Phi(t) = \begin{pmatrix} e^t & 3e^{3t} \\ -e^t & -e^{3t} \end{pmatrix}.\] The method of variation of parameters gives the general solution in the form \[x(t) = \Phi(t)\int \Phi^{-1}(t) g(t) \, dt.\] Direct computations lead to \[\begin{align*} \Phi^{-1}(t) &= \frac{1}{2e^{4t}} \begin{pmatrix} -e^{3t} & -3e^{3t} \\ e^t & e^t \end{pmatrix} = \frac{1}{2}\begin{pmatrix} -e^{-t} & -3e^{-t} \\ e^{-3t} & e^{-3t} \end{pmatrix}. \\ \Phi^{-1}(t)g(t) &= \frac{1}{2}\begin{pmatrix} -e^{-t} & -3e^{-t} \\ e^{-3t} & e^{-3t} \end{pmatrix} \begin{pmatrix}0 \\ t \end{pmatrix} = \begin{pmatrix} -\frac{3}{2}te^{-t} \\ \frac{1}{2}te^{-3t} \end{pmatrix}. \\ \int\Phi^{-1}(t)g(t)\,dt &= \begin{pmatrix} -\frac{3}{2}\int te^{-t} \,dt \\ \frac{1}{2} \int te^{-3t}\, dt \end{pmatrix} = \begin{pmatrix} \frac{3}{2} e^{-t}(t+1) + c_1\\ \frac{1}{18} e^{-3t}(3t+1) + c_2 \end{pmatrix}. \\ \Phi(t)\int\Phi^{-1}(t)g(t)\,dt & = \begin{pmatrix} e^t & 3e^{3t} \\ -e^t & -e^{3t} \end{pmatrix} \begin{pmatrix} \frac{3}{2} e^{-t}(t+1) + c_1\\ \frac{1}{18} e^{-3t}(3t+1) + c_2 \end{pmatrix}. \\ &= \begin{pmatrix} t + \frac{4}{3} \\ -t - \frac{13}{9}\end{pmatrix} + c_1 e^t \begin{pmatrix} 1 \\ -1 \end{pmatrix} + c_2 e^{-3t} \begin{pmatrix} 3 \\ -1 \end{pmatrix} \end{align*}\]

Note the last two terms are the complementary solution. We could have ignored the constants of integration and simply added these on afterwards.