Linear Systems of Ordinary Differential Equations with Constant Coefficients #

Introduction #

Probably the most thoroughly studied set of math models are linear systems of differential equations with constant coefficients (LSDECCs). While most of the world is non-linear, when a linear model is applicable, we should be able to take as great advantage of it as possible. This set of notes walks through the practical aspects of solving linear systems of differential equations with constant coefficients. The emphasis is upon how to solve these problems and in that spirit, proofs are left to references.

The world of homogeneous LSDECC is largely split into two categories: diagonalizable matrices and non-diagonalizable matrices. If a matrix is diagonalizable, then it has real distinct eigenvalues and a distinct eigenvector per eigenvalue. These problems are both easy to identify and simple to solve. Non-diagonalizable matrices are a little tricker to solve, but use the same intuition as the diagonalizable matrices. However, adding a forcing function to the system can quickly make these systems unsolvable. In this work, we will first address diagonalizable matrices. Then the methods will be extended to non-diagonalizable matrices and finally two techniques to handle non-homogeneous LSDECCs will be discussed.

Real Distinct Eigenvalues #

The basic idea is that we want to project the system in to a space that makes it easy to solve, then solve the easy problem and project back into our original space (see Figure 1). In “mathy” language, that means we want to project the system into a new basis (eigenspace) such that the dynamics nicely lay on the axises. So we turn to eigenvalues and eigenvectors. Things get a little tricky when eigenvalues have multiplicity or are complex. In those cases, we create a “subspace” within the eigenbasis such such that we still have orthogonality, but for now we will assume our matrices are nicely diagonalizable.

Consider the equation

where

and

where $a_{ji}\in\mathbb{C}$ . If the systems of equations was a scalar equation, the solution would be an exponential. In the systems of equations, that is also the case:

where $\underline{X}(t)$ is called the fundamental solution matrix. Then the solution to the IVP is,

Now, the question is, what is the exponential operator upon a matrix? There are many different proofs and ways of answering that question that won’t be addressed here. In the case, $\underline{A}$ has real distinct eigenvalues, here’s what to do.

where $\underline{J} = \mathrm{diag}(\lambda_1, \lambda_2,..., \lambda_n)$ is a diagonal matrix of eigenvalues and $\underline{S}$ is is the related eigenbasis (i.e., every eigenvector is a column of $\underline{S}$ ). Thus,

where $e^{\underline{J}t} = \mathrm{diag}(e^{\lambda_1 t}, e^{\lambda_2 t}, ... e^{\lambda_n t})$ . Thus, the solution to the IVP is

Hermitian Matricies #

The derivation above assumes that $\underline{A}$ is diagonalizable. Unfortunately, most systems are not diagonalizable, but it is pragmatic to know when a system is guaranteed to be diagonalizable. The study of diagonalizable systems is actually spectral theory and it often focuses on self-adjoint operators. While that field is too large to cover here, a little bit is quite helpful.

Definition - Hermitian Operator The Hermitian matrix operator ( $\dagger$ ) is the complex conjugate transpose.

Definition - Hermitian Matrix A matrix ( $\underline{H}$ ) is Hermitian if it is equal to it’s Hermitian:

For example,

is Hermitian, but

is not. This leads us to an important theorem.

Theorem 1.1 - A matrix, $A$ , is guaranteed to be diagonalizable (i.e., it has real, distinct eigenvalues) if and only if the matrix is Hermitian.

Corollary 1.1 - The eigenbasis of a Hermitian matrix is a unitary matrix.

Definition - Unitary Matrix A Unitary matrix is a matrix who’s inverse is equal to it’s transpose. That is, $\underline{U}$ is unitary if and only if

The theorem and corollary make identifying and computing the solution to diagonalizable matrices much easier. For reference, a Hermitian matrix is a self-adjoint operator. Self-adjoint operators are important in the Dirac-von Neumann formulation of quantum mechanics where observables are self-adjoint operators and in particular, the Hamiltonian operator is self-adjoint.

Alternative Formulation #

MT note: I think format is more common among physicists and engineers than mathematicians.

Sometimes, authors will move the transform, $\underline{S}$ , around in these equations. So it can get a little tricky and requires some care. One often useful formulation (which only works well for diagonalizable matrices) is to first realize that there are two different projections in the solution equation. First, the exponential solution in the eigenbasis is projected back into the original basis with $\underline{S}e^{\underline{J}t}$ . Second, the initial conditions, $\overline{b}$ , are projected into the eigenbasis with $\underline{S}^{-1}\overline{b}$ . So, now consider setting initial conditions in the eigenbasis as $\overline{c} = \underline{S}^{-1}\overline{b}$ . Then, for the solution may be written as:

where $\overline{u}_n$ is the eigenvector associated with eigenvalue $\lambda_n$ .

Examples #

Example 1 #

Todo

Example 2 #

Use Zero eigen value.

Real Eigenvalues with Multiplicity #

Now suppose we have real eigenvalues, but they are repeated. This means, we do not have distinct eigenvectors per eigenvalue and therefore the matrix is not diagonalizable. But never fear, we can reduce the matrix to a Jordon canonical form. With distinct eigenvalues, the idea is to project into a orthogonal space that decouples the dimensions. Here the idea is the same, except with repeated eigenvalues some dimensions are coupled so the idea is to create a subspace for each of the coupled dimensions such that within the subspace, the projection of the coupled dimensions are orthogonal. Before we discuss the particulars, we must build a little background.

Matrix Exponential #

In the previous discussion, the matrix exponential of the diagonal matrix $\underline{J}$ was given as: $e^{\underline{J}t} = \mathrm{diag}(e^{\lambda_1 t}, e^{\lambda_2 t}, ... e^{\lambda_n t})$ . Here, we would like to give some intuition about where that comes from. Consider the scaler $a$ and the exponential $e^a$ . The exponential can be expanded with a Maclaurin series as

Notice, this definition only depends upon multiplication which is well defined for square matrices. Thus, we can use the Maclaurin series to define the exponential.

Definition - Matrix Exponential Given the matrix $\underline{A} \in \mathbb{C}^{n\times n}$ (see note 1), the exponential of $\underline{A}$ is defined as

Notice that when $\underline{A}$ is diagonal, the series becomes

as given before. However, in most cases, the summation is infinite and nasty, but there are a few properties to help out.

The matrix exponential has the following properties

1. $e^{\underline{A}t}$ is defined $\forall t\in\mathbb{R}$ and $||e^{\underline{A}t}|| \leq ||\underline{A}||\; |t|$ .
2. Given $\underline{A}\underline{B} = \underline{B}\underline{A}$ , $e^{(\underline{A}+\underline{B})} = e^{\underline{A}}e^{\underline{B}} =e^{\underline{B}}e^{\underline{A}}$ .
3. $e^{\underline{A}(t+s)} = e^{\underline{A}t}e^{\underline{A}s} =e^{\underline{A}s}e^{\underline{A}t}\;\forall t,s\in\mathbb{R}$ .
4. $e^{-\underline{A}t}$ exists for all $t\in\mathbb{R}$ .
5. $\frac{d}{dt}e^{\underline{A}t} = \underline{A}e^{\underline{A}t} = e^{\underline{A}t}\underline{A}$ .

While most of these properties are follow from the definition of the scaler exponential, the biggest difference is property (2) which requires $\underline{A}$ and $\underline{B}$ to commute. Forgetting that requirement can easily lead to the wrong answer.

Jordon Canonical Form #

Consider the equation

where $\underline{A} \in \mathbb{R}^{n\times n}$ , $\overline{x}(t): \mathbb{R} \rightarrow \mathbb{R}^n$ and $\overline{b}$ is the initial condition. Suppose $\underline{A}$ has $m$ unique eigenvalues, $\lambda_j$ , such that $\lambda_j$ is repeated $n_j$ times. Thus, $n = \sum_{j=1}^mn_j$ . Repeated eigenvalues means that $n_j$ dimensions are coupled and cannot be separated through projection. We can, however, create an orthogonal subspace to project the $n_j$ dimensions into. First, we must find the generaized eigenvectors. The set of $n_j$ generalized eigenvectors associated with eigenvalue $\lambda_j$ is given by

where $\overline{u}_{j,k=0} = 0$ , $j = 1,2,...,m$ and $k=1,2,...,n_j$ . The trivial case is ignored and $\overline{u}_{j,k}\neq0\;\forall\;k=1,2,...,m$ . Let $\underline{P}_j: n\times n_j$ be the projection from the system basis into the basis formed by generalized eigenvalues associated with eigenvalue $\lambda_j$ defined by

Thus, we can form a generalized eigenbasis with $\underline{P}_j$ as

where $\underline{S}$ is an an $n\times n$ matrix. Using the generalized eigenbasis, the matrix $\underline{A}$ can be composed of

where $\underline{J} = \mathrm{diag}(\underline{B}_1, \underline{B}_2, ...\underline{B}_m)$ is the $n \times n$ Jordon matrix with Jordon blocks $\underline{B}_j: n_j\times n_j$ defined by

where $\underline{I}$ is the $n_j \times n_j$ identity matrix and $\underline{N}$ is the $n_j \times n_j$ nilpotent matrix given by,

Before continuing, let’s clarify the dimensions and indices of everything.

• $n$ - the dimension of $\underline{A}$ and therefore $\underline{J}$ and $\underline{S}$ .
• $m$ - the number of unique eigenvalues of $\underline{A}$ .
• $n_j$ - the multiplicity of the $j^{th}$ eigenvalue of $\underline{A}$ .
• And $n = \sum_{j=1}^mn_j.$

Further notice

• $\underline{P}_j$ is the set of generalized eigenvectors associated with $\lambda_j$ .
• $\underline{B}_j$ is the Jordon block associated with $\lambda_j$ .
• If $n_j=1$ for all $j=1,2,...m$ then the Jordon matrix of the same form as when $\underline{A}$ is diagonalizable.

By applying the decomposition of $\underline{A}$ , the fundamental solution matrix is given by

and the particular solution is

The exponential $e^{\underline{J}t} = \mathrm{diag}(e^{\underline{B}_1t},e^{\underline{B}_2t},...,e^{\underline{B}_mt})$ where

since $\lambda_j\underline{I}$ and $\underline{N}$ commute. Notice that the exponential summation is finite because $\underline{N}$ is nilpotent.

Whew, that is a lot. So what have we done? The matrix $\underline{A}$ can’t be diagonalized so we’ve created subspaces for each set of coupled dimensions using the orthogonality of polynomials which all together form a generalized eigenbasis. We then us the generalized eigenvalues to project into the generalized eigenbasis, solve the system in the generalized eigenbasis, then project back into the original basis.

Complex Eigenvalues - Real Jordon Canonical Form #

If $\underline{A}$ is real, then $\underline{J}$ and $\underline{S}$ may have complex entries, but any complex eigenvalues will occur in conjugate pairs; we can use this fact to simplify multiplicity of complex eigenvalues using the real Jordon canonical form.

Consider the complex eigenvalue pair $\lambda_{ja}, \lambda_{jb}\in\mathbb{C}$ such that $\lambda_{ja} = \lambda_{jb}^* = \alpha_j \pm i\beta_j$ with multiplicity $n_j$ (e.g., for $n_j=3$ , $\lambda_{ja}$ appears three times and $\lambda_{jb}$ appears three times). Rather than using two Jordon Blocks, $\underline{B}_{ja}$ and $\underline{B}_{jb}$ , use one block, $\underline{B}_j^{Re} \in \mathbb{R}^{2n_j \times 2n_j}$ defined as

where $\underline{I}$ is the $2\times2$ identity matrix and

Then the Jordon matrix uses $\underline{B}_j$ for real eigenvalues and $\underline{B}_j^{Re}$ for complex eigenvalues.

Examples #

Example 1 #

Solve the homogeneous systems of differential equations

for initial conditions

Solution From the system of differential equations, we find

Then find the eigenvalues:

Thus, $\lambda = 3$ with multiplicity 2. Next, find generalized eigenvectors

and

Then the generalized eigenbasis is

and the inverse (using Gaussian elimination)

From theory, we know the Jordon matrix should have one block and be

but if we forget (or wish to verify), we can simply use the generalized eigenbasis to find

Then Thus, the fundamental solution matrix is

and the solution to the IVP is

Example 2 #

Multiplicity of three

Example 3 #

Complex eigenvalues and the real form

Non-Homogeneous Linear Systems of ODEs #

Now we can make things more fun by adding in forcing functions. Consider the non-homogenious system of differential equations

where $\underline{A}\in\mathbb{R}^{n\times n}$ , $\overline{x}: \mathbb{R} \rightarrow \mathbb{R}^n$ and $\overline{g}: \mathbb{R} \rightarrow \mathbb{R}^n$ . The function $\overline{g}(t)$ is a forcing function. Just like for scaler ODS’s there are largely two methods for handling forcing functions: method of undetermined coefficients (i.e., guessing) and method of variation of parameters. Unfortunately, the easier, method of undetermined coefficients is far more limited and the difficult (sometimes computationally impossible) method of variation of parameters is the only option.

Method of Undetermined Coefficients #

The method of undetermined coefficients (i.e., guessing) can be quite effective for scalar differential equations, but is rather limited for systems of differential equations. The only guesses that work out nicely are for constants, positive integer powers of $t$ , exponential functions and $\sin$ or $\cos$ .

First, (as always) solve the homogeneous system. This is especially important for exponential functions and $\sin$ or $\cos$ because one or more of the eigenvalues of $\underline{A}$ may conflict with the forcing function. Then, the form of the particular solution can be “guessed” with now vectors for coefficients and back solved into the equation to find the coefficients. Just like for scaler ODEs, the general solution is the sum of the homogeneous and and particular solution

Constants #

The easiest guess is when $\overline{g}(t)$ is constant. The appropriate guess is

where $\overline{a}\in{R}^n$ . Then substitute the particular solution into the non-homogeneous equation to find the constants $\overline{a}$ .

Positive Powers of $t$ #

Positive integer powers of $t$ is also not too tricky. If $\overline{g}(t)$ is order $r$ , then the appropriate guess for the form of the particular solution is,

where $\overline{a}_k\in\mathbb{R}^n$ are the constants to solve by substitution.

Exponential and Sinusoid Functions #

It’s usually easiest to use Euler’s formula to make $\sin$ and $\cos$ into complex exponential functions, so we will consider both exponential and sinusoids together. For an exponential forcing function of the form $\overline{g}(t) = \overline{k}e^{\alpha t}$ , the appropriate guess for the particular solution is

where $\overline{a}\in\mathbb{R}^n$ are the coefficients to back solve, and $r=n_j$ where $n_j$ is the multiplicity of the eigenvalue of the system such that $\lambda_j = \alpha$ . The basic idea is that we must ensure that the particular solution is orthogonal to the eigenspace of the homogeneous solution and we use polynomials to do that. Let’s look at a few cases to understand this better. In the case $\alpha$ is not an eigenvalue of the system (i.e, $n_j=0$ ), then $r=0$ and

In the case $\alpha$ is an eigenvalue of the system with multiplicity $n_j=1$ , then

Notice, by multiplying by $t$ , it ensures that the particular solution is orthogonal to the eigenspace of the homogeneous solution and does not interfere. When $\alpha$ conflicts with eigenvalues of the homogeneous solution of higher order, we simply use a one higher polynomial power as that of the associated Jordon Block to ensure orthogonality.

Method of Variation of Parameters #

Variation of parameters is guaranteed to be valid. However, often the integration is impossible. Let $\underline{A}\in\mathbb{R}^{n\times n}$ and $\overline{g}(t): \mathbb{R} \rightarrow \mathbb{R}^n$ be a continuous map from time into $\mathbb{R}^n$ . Then for every $(t_0, x_0)\in\mathbb{R}^{n+1}$ , the initial value problem

has a unique global solution $x(t,t_0, x_0)$ given by the formula

Examples #

TODO

Notes #

1. Actually, $\underline{A}$ can be any linear operation on a finite normed vector space.

Ref #

[1] Sideris 2013

[2] Kreyszig 2011

[3] Schaeffer 2016

[4] Hirsch 2013