Systems of Linear Differential Equations

When we encounter a situation where more than one dependent variable changes with respect to an independent variable, and these changes are related through a set of differential equations, we call this a system of differential equations. If all the equations in the system are linear, and the dependent variables appear only to the first power with no products of variables, they are called systems of linear differential equations. These systems are fundamental in modeling many real-world phenomena, such as population dynamics, electrical circuits, and mechanical vibrations.

First-Order Linear Systems

The most general form of a first-order linear system of differential equations with $n$ dependent variables $x_1(t), x_2(t), ..., x_n(t)$ and one independent variable $t$ can be written as:

$ \frac{dx_1}{dt} = a_{11}(t)x_1 + a_{12}(t)x_2 + \dots + a_{1n}(t)x_n + f_1(t) $

$ \frac{dx_2}{dt} = a_{21}(t)x_1 + a_{22}(t)x_2 + \dots + a_{2n}(t)x_n + f_2(t) $

...

$ \frac{dx_n}{dt} = a_{n1}(t)x_1 + a_{n2}(t)x_2 + \dots + a_{nn}(t)x_n + f_n(t) $

Here, $a_{ij}(t)$ are known functions of $t$, and $f_i(t)$ are also known functions of $t$. If all $f_i(t)$ are zero, the system is called homogeneous. Otherwise, it is non-homogeneous.

Matrix Representation

This system can be more conveniently expressed using matrix notation. Let $\mathbf{x}(t)$ be the column vector of dependent variables:

$ \mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t) \end{pmatrix} $

Let $\mathbf{A}(t)$ be the coefficient matrix:

$ \mathbf{A}(t) = \begin{pmatrix} a_{11}(t) & a_{12}(t) & \dots & a_{1n}(t) \\ a_{21}(t) & a_{22}(t) & \dots & a_{2n}(t) \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1}(t) & a_{n2}(t) & \dots & a_{nn}(t) \end{pmatrix} $

And let $\mathbf{f}(t)$ be the non-homogeneous term vector:

$ \mathbf{f}(t) = \begin{pmatrix} f_1(t) \\ f_2(t) \\ \vdots \\ f_n(t) \end{pmatrix} $

The derivative of the vector $\mathbf{x}(t)$ is:

$ \frac{d\mathbf{x}}{dt} = \begin{pmatrix} \frac{dx_1}{dt} \\ \frac{dx_2}{dt} \\ \vdots \\ \frac{dx_n}{dt} \end{pmatrix} $

Using these definitions, the system of linear differential equations can be written in matrix form as:

$ \frac{d\mathbf{x}}{dt} = \mathbf{A}(t)\mathbf{x}(t) + \mathbf{f}(t) $

If $\mathbf{f}(t) = \mathbf{0}$, the system is homogeneous:

$ \frac{d\mathbf{x}}{dt} = \mathbf{A}(t)\mathbf{x}(t) $

Solving Homogeneous Systems with Constant Coefficients

A common and important case is when the coefficient matrix $\mathbf{A}$ has constant entries, i.e., $\mathbf{A}(t) = \mathbf{A}$, and the system is homogeneous:

$ \frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x} $

We look for solutions of the form $\mathbf{x}(t) = \mathbf{v}e^{\lambda t}$, where $\mathbf{v}$ is a constant non-zero vector and $\lambda$ is a scalar constant. Substituting this into the differential equation:

$ \frac{d}{dt}(\mathbf{v}e^{\lambda t}) = \mathbf{A}(\mathbf{v}e^{\lambda t}) $

$ \lambda \mathbf{v}e^{\lambda t} = \mathbf{A}\mathbf{v}e^{\lambda t} $

Since $e^{\lambda t}$ is never zero, we can divide by it:

$ \lambda \mathbf{v} = \mathbf{A}\mathbf{v} $

This is an eigenvalue problem. It can be rewritten as:

$ \mathbf{A}\mathbf{v} - \lambda \mathbf{v} = \mathbf{0} $

$ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v} = \mathbf{0} $

where $\mathbf{I}$ is the identity matrix. For a non-trivial solution $\mathbf{v}$ (i.e., $\mathbf{v} \neq \mathbf{0}$), the matrix $(\mathbf{A} - \lambda \mathbf{I})$ must be singular, which means its determinant must be zero:

$ \det(\mathbf{A} - \lambda \mathbf{I}) = 0 $

This equation is called the characteristic equation. The roots $\lambda$ of this equation are the eigenvalues of $\mathbf{A}$, and the corresponding non-zero vectors $\mathbf{v}$ are the eigenvectors.

Case 1: Distinct Real Eigenvalues

If the characteristic equation yields $n$ distinct real eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_n$, each with a corresponding linearly independent eigenvector $\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_n$, then a fundamental set of solutions is given by:

$ \mathbf{x}_1(t) = \mathbf{v}_1 e^{\lambda_1 t}, \quad \mathbf{x}_2(t) = \mathbf{v}_2 e^{\lambda_2 t}, \quad \dots, \quad \mathbf{x}_n(t) = \mathbf{v}_n e^{\lambda_n t} $

The general solution is then a linear combination of these fundamental solutions:

$ \mathbf{x}(t) = c_1 \mathbf{v}_1 e^{\lambda_1 t} + c_2 \mathbf{v}_2 e^{\lambda_2 t} + \dots + c_n \mathbf{v}_n e^{\lambda_n t} $

where $c_1, c_2, \dots, c_n$ are arbitrary constants.

Case 2: Repeated Real Eigenvalues

If the characteristic equation has a real eigenvalue $\lambda$ with algebraic multiplicity $m > 1$, we need to find $m$ linearly independent solutions associated with this eigenvalue.

If for an eigenvalue $\lambda$ with algebraic multiplicity $m$, we can find $m$ linearly independent eigenvectors, then we have $m$ solutions of the form $\mathbf{v}e^{\lambda t}$. However, this is not always the case.

If the geometric multiplicity (the number of linearly independent eigenvectors) is less than the algebraic multiplicity $m$, we need to use generalized eigenvectors. For an eigenvalue $\lambda$ with algebraic multiplicity 2 and geometric multiplicity 1, we find one eigenvector $\mathbf{v}_1$ corresponding to $\lambda$. Then we look for a second linearly independent solution of the form:

$ \mathbf{x}_2(t) = \mathbf{v}_2 t e^{\lambda t} + \mathbf{w} e^{\lambda t} $

where $\mathbf{v}_2$ is a generalized eigenvector and $\mathbf{w}$ is a constant vector satisfying $ (\mathbf{A} - \lambda \mathbf{I})\mathbf{w} = \mathbf{v}_2 $.

A more systematic way to find generalized eigenvectors for an eigenvalue $\lambda$ of algebraic multiplicity $m$ is to solve the system:

$ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v}_1 = \mathbf{0} $ (finds eigenvector $\mathbf{v}_1$)

$ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v}_2 = \mathbf{v}_1 $ (finds generalized eigenvector $\mathbf{v}_2$)

$ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v}_3 = \mathbf{v}_2 $ (finds generalized eigenvector $\mathbf{v}_3$)

...

$ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v}_m = \mathbf{v}_{m-1} $

The $m$ linearly independent solutions associated with $\lambda$ are then:

$ \mathbf{x}_1(t) = \mathbf{v}_1 e^{\lambda t} $

$ \mathbf{x}_2(t) = (\mathbf{v}_1 t + \mathbf{v}_2) e^{\lambda t} $

$ \mathbf{x}_3(t) = \left(\mathbf{v}_1 \frac{t^2}{2!} + \mathbf{v}_2 t + \mathbf{v}_3\right) e^{\lambda t} $

...

$ \mathbf{x}_m(t) = \left(\mathbf{v}_1 \frac{t^{m-1}}{(m-1)!} + \dots + \mathbf{v}_{m-1} t + \mathbf{v}_m\right) e^{\lambda t} $

Case 3: Complex Eigenvalues

If the characteristic equation has a pair of complex conjugate eigenvalues $ \lambda = \alpha \pm i\beta $, where $ \beta \neq 0 $, we can find one complex solution and then extract two real, linearly independent solutions from it.

Let $ \lambda_1 = \alpha + i\beta $ be one complex eigenvalue, and let $ \mathbf{v}_1 = \mathbf{a} + i\mathbf{b} $ be its corresponding eigenvector (where $\mathbf{a}$ and $\mathbf{b}$ are real vectors). The complex solution is:

$ \mathbf{x}_1(t) = (\mathbf{a} + i\mathbf{b}) e^{(\alpha + i\beta)t} $

Using Euler's formula, $ e^{i\theta} = \cos\theta + i\sin\theta $:

$ e^{(\alpha + i\beta)t} = e^{\alpha t} e^{i\beta t} = e^{\alpha t}(\cos(\beta t) + i\sin(\beta t)) $

So, $ \mathbf{x}_1(t) = (\mathbf{a} + i\mathbf{b}) e^{\alpha t}(\cos(\beta t) + i\sin(\beta t)) $

$ \mathbf{x}_1(t) = e^{\alpha t} [ (\mathbf{a}\cos(\beta t) - \mathbf{b}\sin(\beta t)) + i(\mathbf{b}\cos(\beta t) + \mathbf{a}\sin(\beta t)) ] $

The real and imaginary parts of this complex solution are themselves real solutions:

$ \mathbf{x}_{R1}(t) = e^{\alpha t}(\mathbf{a}\cos(\beta t) - \mathbf{b}\sin(\beta t)) $

$ \mathbf{x}_{I1}(t) = e^{\alpha t}(\mathbf{b}\cos(\beta t) + \mathbf{a}\sin(\beta t)) $

These two real solutions are linearly independent. If there is another pair of complex conjugate eigenvalues, we repeat the process.

Solving Non-Homogeneous Systems

For a non-homogeneous system $ \frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x} + \mathbf{f}(t) $, the general solution is the sum of the complementary solution (general solution of the associated homogeneous system) and a particular solution:

$ \mathbf{x}(t) = \mathbf{x}_c(t) + \mathbf{x}_p(t) $

The complementary solution $ \mathbf{x}_c(t) $ is found using the methods described above for homogeneous systems.

Finding a particular solution $ \mathbf{x}_p(t) $ can be done using several methods, including:

1. Method of Undetermined Coefficients

This method is applicable when the non-homogeneous term $ \mathbf{f}(t) $ and its derivatives are of a form that can be represented by a finite set of linearly independent functions (e.g., polynomials, exponentials, sines, cosines, or sums and products of these).

We assume a particular solution $ \mathbf{x}_p(t) $ that has the same form as $ \mathbf{f}(t) $, with undetermined coefficients. We then substitute this assumed solution into the non-homogeneous differential equation and solve for the coefficients.

If any term in the assumed $ \mathbf{x}_p(t) $ is a solution to the homogeneous system, we must modify the assumed form by multiplying by $ t $ (or $ t^k $ if necessary, where $ k $ is the smallest integer such that no term in the modified form is a homogeneous solution).

2. Variation of Parameters

This method is more general and can be used for any $ \mathbf{f}(t) $. It is particularly useful when the non-homogeneous term is complex or when the coefficients of the matrix $\mathbf{A}$ are not constant.

Consider the homogeneous system $ \frac{d\mathbf{x}}{dt} = \mathbf{A}(t)\mathbf{x} $. Let $ \mathbf{x}_1(t), \dots, \mathbf{x}_n(t) $ be a fundamental set of solutions. We can form a fundamental matrix $ \boldsymbol{\Phi}(t) $ whose columns are these solutions: $ \boldsymbol{\Phi}(t) = [\mathbf{x}_1(t) | \mathbf{x}_2(t) | \dots | \mathbf{x}_n(t)] $. It satisfies $ \boldsymbol{\Phi}'(t) = \mathbf{A}(t)\boldsymbol{\Phi}(t) $.

For the non-homogeneous system $ \frac{d\mathbf{x}}{dt} = \mathbf{A}(t)\mathbf{x} + \mathbf{f}(t) $, we assume a particular solution of the form:

$ \mathbf{x}_p(t) = \boldsymbol{\Phi}(t)\mathbf{u}(t) $

where $ \mathbf{u}(t) $ is a vector of unknown functions $ \mathbf{u}(t) = \begin{pmatrix} u_1(t) \\ \vdots \\ u_n(t) \end{pmatrix} $.

Taking the derivative: $ \mathbf{x}_p'(t) = \boldsymbol{\Phi}'(t)\mathbf{u}(t) + \boldsymbol{\Phi}(t)\mathbf{u}'(t) $.

Substitute into the non-homogeneous equation:

$ \boldsymbol{\Phi}'(t)\mathbf{u}(t) + \boldsymbol{\Phi}(t)\mathbf{u}'(t) = \mathbf{A}(t)\boldsymbol{\Phi}(t)\mathbf{u}(t) + \mathbf{f}(t) $

Since $ \boldsymbol{\Phi}'(t) = \mathbf{A}(t)\boldsymbol{\Phi}(t) $, we have:

$ \mathbf{A}(t)\boldsymbol{\Phi}(t)\mathbf{u}(t) + \boldsymbol{\Phi}(t)\mathbf{u}'(t) = \mathbf{A}(t)\boldsymbol{\Phi}(t)\mathbf{u}(t) + \mathbf{f}(t) $

This simplifies to:

$ \boldsymbol{\Phi}(t)\mathbf{u}'(t) = \mathbf{f}(t) $

To find $ \mathbf{u}'(t) $, we multiply by the inverse of the fundamental matrix, $ \boldsymbol{\Phi}^{-1}(t) $:

$ \mathbf{u}'(t) = \boldsymbol{\Phi}^{-1}(t)\mathbf{f}(t) $

Now, we integrate to find $ \mathbf{u}(t) $:

$ \mathbf{u}(t) = \int \boldsymbol{\Phi}^{-1}(t)\mathbf{f}(t) dt $

Finally, the particular solution is $ \mathbf{x}_p(t) = \boldsymbol{\Phi}(t)\mathbf{u}(t) $.

Systems of Higher-Order Linear Differential Equations

A system of higher-order linear differential equations can be converted into an equivalent system of first-order linear differential equations. This is done by introducing new dependent variables for the derivatives of the original variables.

Consider a system involving two variables $x(t)$ and $y(t)$:

$ a_{11} \frac{d^2x}{dt^2} + a_{12} \frac{dy}{dt} + a_{13} x = f_1(t) $ $ a_{21} \frac{dx}{dt} + a_{22} \frac{d^2y}{dt^2} + a_{23} y = f_2(t) $

We introduce new variables:

$ x_1 = x, \quad x_2 = \frac{dx}{dt} $ $ y_1 = y, \quad y_2 = \frac{dy}{dt} $

Then the system becomes a system of first-order equations:

$ \frac{dx_1}{dt} = x_2 $ $ a_{11} \frac{dx_2}{dt} = -a_{12} y_2 - a_{13} x_1 + f_1(t) \implies \frac{dx_2}{dt} = -\frac{a_{12}}{a_{11}} y_2 - \frac{a_{13}}{a_{11}} x_1 + \frac{f_1(t)}{a_{11}} $ $ \frac{dy_1}{dt} = y_2 $ $ a_{22} \frac{dy_2}{dt} = -a_{21} x_2 - a_{23} y_1 + f_2(t) \implies \frac{dy_2}{dt} = -\frac{a_{21}}{a_{22}} x_2 - \frac{a_{23}}{a_{22}} y_1 + \frac{f_2(t)}{a_{22}} $

This new system can then be written in matrix form and solved using the methods for first-order systems. This reduction is crucial for applying the powerful tools of linear algebra to higher-order systems.

Applications of Systems of Linear Differential Equations

Systems of linear differential equations are used to model a wide array of phenomena:

• Electrical Circuits: Analyzing circuits with multiple loops and components (resistors, inductors, capacitors) leads to systems of differential equations describing the voltage and current in different parts of the circuit.
• Mechanical Vibrations: Coupled oscillators, such as two masses connected by springs, exhibit behavior described by systems of differential equations. The motion of each mass influences the others.
• Population Dynamics: Predator-prey models (e.g., Lotka-Volterra equations, though often non-linear) and competition models between species can be approximated by linear systems, especially for small populations or when analyzing deviations from an equilibrium.
• Chemical Reactions: Modeling the rates of change of concentrations of multiple reactants and products in a chemical process.
• Control Systems: Designing controllers for complex systems often involves analyzing the stability and response of linear systems.