Linear systems

First-order systems, matrix exponentials and transform methods.

In this chapter 4

Tick a unit to add it to your PDF basket.

Systems of linear equations

Systems of linear equations

Method

Linear systems and the matrix exponential

Adapted

Add to PDF

Gives the solution of the first-order linear system x' = A(t)x in terms of a fundamental matrix, shows that for constant A the fundamental matrix is exp(tA) and the unique solution of the initial-value problem is exp(tA)x_0, and turns this into a recipe that reads the solution off the eigenvalues and eigenvectors of A, including the two real solutions carried by each complex conjugate pair.

Prerequisites: Functions of matrices and the matrix exponential, Structure of linear equations and the Wronskian, Eigenvalues and eigenvectors, Jordan normal form, Existence and uniqueness of solutions, Determinants, Linear systems and rank, Generalized eigenproblems and normal modes

Definition(The system, its fundamental matrices and its propagator)AdaptedAdapted from Dummit & Foote (2004), §12.3, the discussion between Exercises 52 and 53, pp. 505-506

Fix an interval $I \subseteq \mathbb{R}$ together with an $n \times n$ matrix $A(t)$ whose entries are continuous complex-valued functions of $t \in I$ . The homogeneous linear system is

x'(t) = A(t)\, x(t) ,

where a solution is a differentiable function $x : I \to \mathbb{C}^n$ satisfying this identity at every $t \in I$ ; derivatives and integrals of vector- and matrix-valued functions are taken entry by entry. When $A$ does not depend on $t$ the system is said to have constant coefficients, and we write it $x' = Ax$ .

An $n \times n$ matrix-valued function $\Phi$ on $I$ whose $n$ columns are linearly independent solutions of the system is called a fundamental matrix of it. Equivalently — reading the defining equation column by column — $\Phi$ satisfies the matrix differential equation

\Phi'(t) = A(t)\, \Phi(t)

and its columns are linearly independent as functions on $I$ . If moreover $\Phi(t_0) = I_n$ for a chosen $t_0 \in I$ , we call $\Phi$ the principal fundamental matrix at $t_0$ . (Throughout the unit $I_m$ denotes the $m \times m$ identity matrix, the letter $I$ without a subscript being reserved for the interval.) The two-variable matrix

U(t, t_0) = \Phi(t)\, \Phi(t_0)^{-1}

built from any fundamental matrix is called the propagator (or state-transition matrix) of the system: the statements below show that it does not depend on which $\Phi$ was used, and that it carries the state at time $t_0$ to the state at time $t$ .

Method(Existence, uniqueness and the dimension of the solution space)AdaptedAdapted from Coddington (1989), Ch. 6, §7, Theorem 6, p. 256; Dummit & Foote (2004), §12.3, the discussion preceding Exercise 53, p. 506

Let $A(t) = [a_{jk}(t)]$ and $b(t) = (b_1(t), \dots, b_n(t))$ have entries continuous on an interval $I$ containing $t_0$ . For every $x_0 \in \mathbb{C}^n$ the initial-value problem

x' = A(t)x + b(t), \qquad x(t_0) = x_0

has one and only one solution on the whole of $I$ ((Coddington (1989), Ch. 6, §7, Theorem 6, p. 256)). Three consequences, for the homogeneous system $b \equiv 0$ :

The solution space. The set of solutions of $x' = A(t)x$ on $I$ is a vector space of dimension exactly $n$ over $\mathbb{C}$ , and for each fixed $t_1 \in I$ the evaluation map $x \mapsto x(t_1)$ is a linear isomorphism of it onto $\mathbb{C}^n$ .
Fundamental matrices are nonsingular everywhere or nowhere. If $\Phi$ solves $\Phi' = A(t)\Phi$ on $I$ , then either $\det \Phi(t) \neq 0$ for every $t \in I$ , in which case $\Phi$ is a fundamental matrix, or $\det \Phi(t) = 0$ for every $t \in I$ . In particular it is enough to test one point.
The general solution and the propagator. If $\Phi$ is a fundamental matrix then every solution is $x(t) = \Phi(t)c$ for exactly one constant vector $c \in \mathbb{C}^n$ , the solution with $x(t_0) = x_0$ being $x(t) = \Phi(t)\Phi(t_0)^{-1}x_0$ . Consequently the propagator $U(t,t_0) = \Phi(t)\Phi(t_0)^{-1}$ is the same for every fundamental matrix $\Phi$ , and $\Phi Q$ is again a fundamental matrix for every nonsingular constant matrix $Q$ ((Dummit & Foote (2004), §12.3, Exercise 54, p. 506)).

The interval is not shrunk anywhere: linear systems have non-local existence, and the solution lives as long as the coefficients are continuous. Coddington reaches it by observing that the successive-approximation proof of his local existence theorem, which is stated for an interval $|t - t_0| \le a$ , goes through unchanged on an arbitrary interval $\alpha \le t \le \beta$ containing $t_0$ ((Coddington (1989), Ch. 6, §7, p. 256)); $I$ may even be infinitely long.

Idea of the proofOriginal

Everything here is the single sentence "the solution is determined by its value at one instant", which is what existence and uniqueness say, read in three directions. Read as a statement about a linear map it says that evaluation at $t_1$ is injective (uniqueness) and surjective (existence) from the solution space onto $\mathbb{C}^n$ , which forces the dimension to be $n$ . Read as a statement about a matrix of solutions it says that a nontrivial linear relation between the columns at one instant propagates to every instant, because the solution that the relation builds vanishes at that instant and therefore vanishes identically. Read as a formula it says that the only freedom in a solution is the constant vector $c$ of coefficients, and matching $\Phi(t_0)c$ to $x_0$ fixes $c$ .

ProofOriginal

We take the existence-and-uniqueness theorem as given and derive parts 1-3 from it.

The solution space. The set $S$ of solutions on $I$ is a complex vector space: if $x$ and $y$ solve the system and $\alpha, \beta \in \mathbb{C}$ , then $(\alpha x + \beta y)' = \alpha A x + \beta A y = A(\alpha x + \beta y)$ , because differentiation is linear and $A(t)$ acts linearly. Fix $t_1 \in I$ and let $E : S \to \mathbb{C}^n$ be $E(x) = x(t_1)$ , which is linear. $E$ is surjective: given $v \in \mathbb{C}^n$ , the theorem applied with $t_0 = t_1$ , $x_0 = v$ produces a solution on all of $I$ with that value. $E$ is injective: if $x \in S$ and $x(t_1) = 0$ , then $x$ and the zero function are both solutions of the initial-value problem with data $(t_1, 0)$ — the zero function solves it because $0' = 0 = A(t)\cdot 0$ — so uniqueness gives $x \equiv 0$ . A linear isomorphism onto $\mathbb{C}^n$ forces $\dim S = n$ .
Nonsingular everywhere or nowhere. Suppose $\det \Phi(t_1) = 0$ for some $t_1 \in I$ . Then the columns of $\Phi(t_1)$ are linearly dependent, so there is $c \neq 0$ with $\Phi(t_1)c = 0$ . The function $y(t) = \Phi(t)c$ is a solution of $y' = A(t)y$ , since $y' = \Phi' c = A\Phi c = Ay$ — here $c$ being constant is what lets it pass through the derivative — and $y(t_1) = 0$ . By the injectivity just proved, $y \equiv 0$ , that is $\Phi(t)c = 0$ for every $t$ with the same $c \neq 0$ . Hence $\det \Phi(t) = 0$ for every $t \in I$ . Contrapositively, if $\det$ is nonzero at one point it is nonzero at all of them; and then at each $t$ the columns of $\Phi(t)$ are independent vectors, so the columns of $\Phi$ are independent as functions and $\Phi$ is a fundamental matrix. (Conversely, if the columns are independent as functions, they cannot be dependent at any single point, by the argument just given, so the determinant never vanishes. The two descriptions of a fundamental matrix in the definition therefore agree.)
The general solution. Let $\Phi$ be a fundamental matrix and $x$ any solution. Put $c = \Phi(t_1)^{-1}x(t_1)$ for some fixed $t_1 \in I$ , legitimate by part 2. Then $\Phi(t)c$ and $x(t)$ are solutions agreeing at $t_1$ , so they agree on $I$ by injectivity again; and $c$ is unique because $\Phi(t_1)c = \Phi(t_1)\tilde c$ with $\Phi(t_1)$ invertible gives $c = \tilde c$ . Choosing $t_1 = t_0$ and $c = \Phi(t_0)^{-1}x_0$ gives the stated formula. If $\Psi$ is another fundamental matrix, each of its columns is $\Phi(\cdot)c_j$ for some $c_j$ , so $\Psi = \Phi Q$ with $Q$ the constant matrix of columns $c_j$ , and $Q = \Phi(t_0)^{-1}\Psi(t_0)$ is invertible; therefore

\Psi(t)\Psi(t_0)^{-1} = \Phi(t) Q Q^{-1} \Phi(t_0)^{-1} = \Phi(t)\Phi(t_0)^{-1} ,

so the propagator is independent of the choice. Running this backwards, for any nonsingular constant $Q$ the matrix $\Phi Q$ satisfies $(\Phi Q)' = \Phi' Q = A \Phi Q$ and has determinant $\det \Phi \cdot \det Q \neq 0$ , so it is a fundamental matrix. $\square$

Method(Constant coefficients: the exponential is the propagator)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.4, the Theorem on p. 129; Dummit & Foote (2004), §12.3, Exercise 53, p. 506

Let $A$ be a constant $n \times n$ complex matrix. Then $t \mapsto \exp(tA)$ is a fundamental matrix of $x' = Ax$ on all of $\mathbb{R}$ ; it is the principal one at $t = 0$ . The unique solution of

x' = Ax , \qquad x(0) = x_0

is $x(t) = \exp(tA)\,x_0$ ((Hirsch, Smale & Devaney (2013), §6.4, Theorem, p. 129)), and the propagator depends only on the elapsed time:

U(t, t_0) = \exp\bigl((t - t_0)A\bigr) .

Three readings of the same fact are worth separating. (i) Algebraically, $\exp(tA)$ is invertible with inverse $\exp(-tA)$ , so the state can be propagated backwards as well as forwards. (ii) Dynamically, $U(t+s, t_0+s) = U(t, t_0)$ : the system is autonomous, and only elapsed time matters. (iii) Structurally, $\exp\bigl((t+s)A\bigr) = \exp(tA)\exp(sA)$ , so $\{\exp(tA)\}_{t \in \mathbb{R}}$ is a one-parameter group of matrices whose generator is $A$ ; this is the finite-dimensional case of the semigroup picture of evolution equations.

Idea of the proofOriginal

All of it comes from the one differentiation rule $\frac{d}{dt}\exp(tA) = A\exp(tA)$ established in the unit on the matrix exponential: that rule is literally the statement $\Phi' = A\Phi$ for $\Phi(t) = \exp(tA)$ , and $\exp(0 \cdot A) = I_n$ is nonsingular, so the previous statement does the rest with no new work. The time-translation form of the propagator then follows from the product rule for exponentials of commuting matrices, since $tA$ and $-t_0A$ always commute.

ProofAdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.4, the Theorem on p. 129

$\exp(tA)$ is a fundamental matrix. By the properties of the exponential (matrix exponential) the matrix-valued function $\Phi(t) = \exp(tA)$ is differentiable with $\Phi'(t) = A\exp(tA) = A\Phi(t)$ , so each of its columns solves $x' = Ax$ . At $t = 0$ the series collapses to its first term, $\Phi(0) = I_n$ , which is nonsingular, so by part 2 of the previous statement $\Phi$ is a fundamental matrix, and by construction it is principal at $t = 0$ .
The initial-value problem. By part 3 of the previous statement with $t_0 = 0$ , the unique solution with $x(0) = x_0$ is $\Phi(t)\Phi(0)^{-1}x_0 = \exp(tA)x_0$ . This is the theorem of (Hirsch, Smale & Devaney (2013), §6.4, p. 129), whose own proof is self-contained: it differentiates $\exp(tA)x_0$ to verify that it is a solution, and for uniqueness sets $z(t) = \exp(-tA)y(t)$ for a second solution $y$ and computes $z' \equiv 0$ , so that $z$ is the constant $x_0$ . That route needs only the differentiation rule and the invertibility of the exponential, not the dimension count of the previous statement.
The propagator. The matrices $tA$ and $-t_0A$ commute, so the product rule for exponentials of commuting matrices gives

U(t,t_0) = \exp(tA)\exp(t_0A)^{-1} = \exp(tA)\exp(-t_0A) = \exp\bigl((t-t_0)A\bigr) ,

using also $\exp(t_0A)^{-1} = \exp(-t_0A)$ . Replacing $(t, t_0)$ by $(t+s, t_0+s)$ changes neither $t - t_0$ nor therefore $U$ ; and the same identity with $t_0$ replaced by $-s$ reads $\exp\bigl((t+s)A\bigr) = \exp(tA)\exp(sA)$ . $\square$

Method(Reading the solution off the spectrum)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.1, p. 108 and the Theorem on p. 111

Let $A$ be a constant $n \times n$ matrix.

Eigensolutions. If $Av = \lambda v$ with $v \neq 0$ , then $x(t) = e^{\lambda t}v$ solves $x' = Ax$ . Conversely every solution of the form $e^{\lambda t}v$ with $v \neq 0$ arises this way.
Diagonalizable $A$ . If $A$ has eigenvectors $v_1, \dots, v_n$ forming a basis of $\mathbb{C}^n$ , with $Av_j = \lambda_j v_j$ , then

x(t) = \sum_{j=1}^{n} c_j e^{\lambda_j t} v_j , \qquad c \in \mathbb{C}^n ,

is the general solution, and $c = S^{-1}x_0$ where $S$ is the matrix with columns $v_j$ . This covers in particular the case of $n$ distinct eigenvalues, in which Hirsch, Smale and Devaney obtain exactly this formula by diagonalizing and solving the decoupled system ((Hirsch, Smale & Devaney (2013), §6.1, p. 108)). 3. Real form for a complex pair. Let $A$ be real, let $\lambda = \alpha + i\beta$ with $\beta \neq 0$ be an eigenvalue and $v = p + iq$ (with $p, q \in \mathbb{R}^n$ ) a corresponding eigenvector. Then

x_{\mathrm{re}}(t) = e^{\alpha t}\bigl( p\cos\beta t - q \sin\beta t \bigr) , \qquad x_{\mathrm{im}}(t) = e^{\alpha t}\bigl( p\sin\beta t + q\cos\beta t\bigr)

are real solutions of $x' = Ax$ , they are linearly independent over $\mathbb{R}$ — and hence, both being real-valued, over $\mathbb{C}$ as well — and their complex span is the same as the complex span of the two complex eigensolutions belonging to $\lambda$ and $\bar\lambda$ . They are the real and imaginary parts of $e^{\lambda t}v$ , and this is the device by which real solutions are extracted from a complex one in (Hirsch, Smale & Devaney (2013), §3.2, pp. 44-47). 4. Defective $A$ . If $A$ is not diagonalizable, parts 1-3 do not supply enough solutions, and the missing ones are obtained from $\exp(tA)$ by the splitting $A = D + N$ with $D$ diagonalizable, $N$ nilpotent and $DN = ND$ (matrix exponential). Every entry of $\exp(tA)$ is then a sum of terms $(\text{polynomial in } t) \times e^{\lambda t}$ over the eigenvalues $\lambda$ , the degree of the polynomial attached to $\lambda$ being at most one less than the size of the largest Jordan block for $\lambda$ .

Hirsch, Smale and Devaney state part 2 together with part 3 in a single theorem for the case of $n$ distinct eigenvalues, $k_1$ of them real and $k_2$ conjugate pairs: the change of coordinates $T$ built from the eigenvectors (and from the real and imaginary parts of the complex ones) puts $A$ into a block form with the real eigenvalues on the diagonal and $2\times 2$ blocks $\begin{pmatrix} \alpha_j & \beta_j \\ -\beta_j & \alpha_j\end{pmatrix}$ , and the general solution is $T$ applied to the explicit column of exponentials, sines and cosines displayed there ((Hirsch, Smale & Devaney (2013), §6.1, Theorem, p. 111)). The version above separates the two hypotheses, because a diagonalizable matrix with repeated eigenvalues — a real symmetric matrix, for instance — is covered by part 2 but not by a theorem that assumes the eigenvalues distinct.

Idea of the proofOriginal

A single eigendirection reduces the system to one scalar equation, because along $v$ the matrix $A$ only multiplies by $\lambda$ ; so the whole of parts 1 and 2 is the scalar fact $x' = \lambda x \Rightarrow x = ce^{\lambda t}$ , applied once per eigendirection and reassembled. Part 3 is bookkeeping with conjugation: a real matrix maps real vectors to real vectors, so it commutes with taking real and imaginary parts, and therefore the real and imaginary parts of any complex solution are separately solutions. The two real solutions so obtained are not new information — they span the same thing as the conjugate pair of complex eigensolutions — but they are the form in which an answer to a real physical problem has to be delivered.

ProofOriginal

Eigensolutions. If $Av = \lambda v$ , then $x(t) = e^{\lambda t}v$ has $x'(t) = \lambda e^{\lambda t} v = e^{\lambda t}(\lambda v) = e^{\lambda t}Av = A x(t)$ , using the linearity of $A$ to move the scalar $e^{\lambda t}$ across it. Conversely if $x(t) = e^{\lambda t}v$ with $v \neq 0$ solves the system, then differentiating and cancelling the nonzero scalar $e^{\lambda t}$ gives $\lambda v = Av$ , so $\lambda$ is an eigenvalue with eigenvector $v$ .
Diagonalizable $A$ . Each $e^{\lambda_j t}v_j$ is a solution by part 1, so by linearity every $\sum_j c_j e^{\lambda_j t}v_j$ is one. These exhaust the solutions by a dimension count: let $\Phi(t)$ be the matrix with columns $e^{\lambda_j t}v_j$ . Then $\Phi(0) = S$ , which is nonsingular because the $v_j$ form a basis, so $\Phi$ is a fundamental matrix by part 2 of the existence statement, and by part 3 of it every solution is $\Phi(t)c$ — which is the displayed sum. Setting $t = 0$ gives $x_0 = Sc$ , hence $c = S^{-1}x_0$ .
Real form for a complex pair. Expanding $e^{\lambda t}v$ with Euler's formula and $v = p + iq$ ,

e^{(\alpha + i\beta)t}(p + iq) = e^{\alpha t}\bigl(\cos\beta t + i \sin\beta t\bigr)(p + iq) = x_{\mathrm{re}}(t) + i\,x_{\mathrm{im}}(t)

with $x_{\mathrm{re}}, x_{\mathrm{im}}$ exactly as displayed, and both are real because $p$ , $q$ , $\cos \beta t$ and $\sin\beta t$ are. Write $z(t) = e^{\lambda t}v$ , a solution by part 1. Since $A$ has real entries, $A\overline{w} = \overline{Aw}$ for every $w \in \mathbb{C}^n$ — this is where the hypothesis that $A$ is real is used, and it is the only place — so taking complex conjugates in $z' = Az$ gives $\overline{z}' = A\overline{z}$ : the conjugate is a solution too. Hence so are the combinations

x_{\mathrm{re}} = \tfrac12(z + \overline z), \qquad x_{\mathrm{im}} = \tfrac1{2i}(z - \overline z) ,

and these are the real and imaginary parts. The same two formulas, read backwards as $z = x_{\mathrm{re}} + i x_{\mathrm{im}}$ and $\overline z = x_{\mathrm{re}} - ix_{\mathrm{im}}$ , show that the complex span of $\{x_{\mathrm{re}}, x_{\mathrm{im}}\}$ equals that of $\{z, \overline z\}$ ; and $\overline{z}(t) = e^{\bar\lambda t}\overline v$ is the eigensolution for $\bar\lambda$ , since conjugating $Av = \lambda v$ gives $A \overline v = \bar\lambda \overline v$ .

For linear independence over $\mathbb{R}$ , suppose $a\,x_{\mathrm{re}} + b\,x_{\mathrm{im}} = 0$ on an interval with $a, b$ real. At $t = 0$ this reads $ap + bq = 0$ . Differentiating and evaluating at $t = 0$ , using $x_{\mathrm{re}}'(0) = \alpha p - \beta q$ and $x_{\mathrm{im}}'(0) = \alpha q + \beta p$ , gives $a(\alpha p - \beta q) + b(\alpha q + \beta p) = 0$ ; subtracting $\alpha$ times the first relation leaves $\beta(bp - aq) = 0$ , and $\beta \neq 0$ gives $bp - aq = 0$ . Combining,

a(ap + bq) + b(bp - aq) = (a^2 + b^2)\,p = 0, \qquad b(ap+bq) - a(bp - aq) = (a^2+b^2)\,q = 0 .

Now $p$ and $q$ cannot both be zero, since $v = p + iq \neq 0$ ; so $a^2 + b^2 = 0$ and hence $a = b = 0$ .

Independence over $\mathbb{C}$ follows from the real case because both functions take values in $\mathbb{R}^n$ : if $a\,x_{\mathrm{re}} + b\,x_{\mathrm{im}} = 0$ with $a, b \in \mathbb{C}$ , splitting that relation into its real and imaginary parts gives the two real relations $(\operatorname{Re}a)\,x_{\mathrm{re}} + (\operatorname{Re}b)\,x_{\mathrm{im}} = 0$ and $(\operatorname{Im}a)\,x_{\mathrm{re}} + (\operatorname{Im}b)\,x_{\mathrm{im}} = 0$ , whose coefficients all vanish by the paragraph above; hence $a = b = 0$ . This is what makes the comparison of complex spans in the previous paragraph a comparison of two-dimensional spaces.

Defective $A$ . This is the content of the method for computing $\exp(tA)$ in the matrix exponential, combined with the statement above that $\exp(tA)$ is a fundamental matrix; nothing further is proved here. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Continuity of $A$ on the whole interval is what makes existence non-local. For the scalar system $x' = \frac{1}{1-t}x$ on $I = (-\infty, 1)$ the solution through $x(0) = 1$ is $x(t) = (1-t)^{-1}$ , which blows up as $t \to 1^-$ : the interval on which the coefficient is continuous is exactly the interval of existence, and no larger one is claimed.
Linearity is what rules out blow-up inside $I$ and what makes the solution set a vector space. The nonlinear scalar equation $x' = x^2$ with $x(0) = 1$ has solution $(1-t)^{-1}$ although its right-hand side is continuous everywhere; so the non-local conclusion of the existence-and-uniqueness theorem — that the solution lives on the whole interval on which the coefficients are continuous — is not a consequence of continuity alone, and linearity in $x$ is what buys it.
$A$ real is used only in part 3 of the spectral statement, to conjugate the equation. For complex $A$ an eigenvalue need not come with its conjugate, and there is no real form to extract.
$\beta \neq 0$ in part 3 is needed for the independence of $x_{\mathrm{re}}$ and $x_{\mathrm{im}}$ : at $\beta = 0$ they become $e^{\alpha t}p$ and $e^{\alpha t}q$ , which are dependent whenever $p$ and $q$ are parallel — as they may be for a real eigenvalue, where $v$ can be scaled to be real and $q = 0$ .
A basis of eigenvectors is what part 2 needs, not distinct eigenvalues. Conversely, distinctness is sufficient but is also strictly stronger: $A = I_2$ has a repeated eigenvalue and is diagonalizable, while $A = \begin{pmatrix} 0&1\\0&0\end{pmatrix}$ has a repeated eigenvalue and is not, and for it $\exp(tA) = \begin{pmatrix}1&t\\0&1\end{pmatrix}$ carries a solution $(t, 1)^{\mathsf T}$ that is not of the form $e^{\lambda t}v$ at all. That is why part 4 cannot be dispensed with.
$\Phi$ nonsingular at one point is all that part 2 of the existence statement requires, and the conclusion genuinely uses that $\Phi$ solves the differential equation: an arbitrary matrix-valued function can of course be singular at some points and not others. For $n = 1$ this is the familiar fact that a nontrivial solution of a linear scalar equation has no zeros.

Method(Solving x' = Ax with constant A)Original

Compute the characteristic polynomial and its roots. $p(\lambda) = \det(\lambda I_n - A)$ ; its roots are the eigenvalues, with algebraic multiplicities (eigenvalues and eigenvectors). For $2\times2$ matrices use $p(\lambda) = \lambda^2 - (\operatorname{tr}A)\lambda + \det A$ .
For each eigenvalue $\lambda$ , solve $(A - \lambda I_n)v = 0$ and record the dimension of the kernel (the geometric multiplicity). If the geometric multiplicities add up to $n$ , $A$ is diagonalizable and steps 3-4 finish the job; otherwise go to step 5.
Real eigenvalue $\lambda$ : each basis vector $v$ of its eigenspace contributes the solution $e^{\lambda t}v$ .
Complex pair $\lambda = \alpha \pm i\beta$ , $\beta > 0$ , with $A$ real: take an eigenvector $v = p + iq$ for $\alpha + i\beta$ only — the conjugate eigenvalue contributes nothing new — and write down the two real solutions

e^{\alpha t}(p\cos\beta t - q\sin\beta t), \qquad e^{\alpha t}(p\sin\beta t + q\cos\beta t) .

Do not also use $\alpha - i\beta$ : that would double-count. 5. Defective eigenvalue: either compute $\exp(tA)$ by the splitting $A = D + N$ of the matrix exponential and take its columns, or build a chain of generalized eigenvectors $v_1, \dots, v_m$ with $(A - \lambda I_n)v_1 = 0$ and $(A-\lambda I_n)v_{k+1} = v_k$ , which contributes the $m$ solutions

e^{\lambda t}\sum_{j=0}^{k-1} \frac{t^{j}}{j!}\, v_{k-j} \qquad (k = 1, \dots, m)

— that is, $e^{\lambda t}v_1$ , then $e^{\lambda t}(v_2 + t v_1)$ , then $e^{\lambda t}(v_3 + tv_2 + \tfrac12 t^2 v_1)$ , and so on. 6. Assemble the fundamental matrix $\Phi(t)$ from the $n$ solutions as columns, and check $\det\Phi(t) \neq 0$ at the single convenient point $t = 0$ . 7. Fit the initial condition by solving the linear system $\Phi(0)c = x_0$ (linear systems); the answer is $x(t) = \Phi(t)c$ . Equivalently $x(t) = \Phi(t)\Phi(0)^{-1}x_0$ , and $\exp(tA) = \Phi(t)\Phi(0)^{-1}$ if the exponential itself is wanted. 8. Two cheap checks. $x(0)$ must equal $x_0$ , and $x'(0)$ must equal $Ax_0$ ; and $\det \exp(tA) = e^{t \operatorname{tr}A}$ , so $\det\Phi(t) = \det\Phi(0)\,e^{t\operatorname{tr}A}$ (determinants).

Reading the answer without computing it. The real parts of the eigenvalues decide the fate of every solution: $\operatorname{Re}\lambda < 0$ for all $\lambda$ means every solution decays to $0$ ; one eigenvalue with $\operatorname{Re}\lambda > 0$ means almost every solution grows; $\operatorname{Im}\lambda = \pm\beta$ contributes oscillation at angular frequency $\beta$ , that is period $2\pi/\beta$ . Eigenvalues of very different magnitude make the system stiff and its numerical solution delicate (stiffness); the planar case is classified completely in the classification of planar linear systems.

Remark(Scalar equations are systems)AdaptedAdapted from Coddington (1989), Ch. 6, §8, pp. 258-259; Dummit & Foote (2004), §12.3, Exercise 57 and the discussion preceding it, p. 507

A scalar equation of order $n$ ,

y^{(n)} + a_1 y^{(n-1)} + \dots + a_{n-1}y' + a_n y = 0 ,

becomes a first-order system in the vector of derivatives: setting $x = (y, y', \dots, y^{(n-1)})^{\mathsf T}$ , the first $n-1$ components of $x' = Ax$ restate the definitions $x_k' = x_{k+1}$ , and the last restates the equation. Coddington uses exactly this substitution to transfer existence and uniqueness from systems to $n$ -th order equations ((Coddington (1989), Ch. 6, §8, Theorems 7 and 8, pp. 259-260)); Dummit and Foote note that the resulting coefficient matrix is the transpose of the companion matrix of the characteristic polynomial, and that the first row of a fundamental matrix of the system is a basis of solutions of the scalar equation ((Dummit & Foote (2004), §12.3, Exercise 57 and the preceding discussion, p. 507)).

Two consequences are worth keeping. The characteristic polynomial of the companion matrix is the characteristic polynomial of the scalar equation, so "eigenvalues of $A$ " and "roots of the indicial polynomial" of the constant-coefficient equation are the same list, and a root of multiplicity $m$ gives a single Jordan block of size $m$ . (That last point is ours, not the sources': for the matrix $A$ above, deleting the first column and the last row of $A - \lambda I_n$ leaves a triangular matrix with $1$ 's on the diagonal, so $\operatorname{rank}(A - \lambda I_n) \ge n-1$ and every eigenspace is one-dimensional — a companion matrix is maximally defective, and step 5 of the method is unavoidable for a scalar equation with a repeated root.) And the Wronskian of the theory of linear equations is the determinant of the corresponding fundamental matrix, which is why the two nonvanishing criteria match. In the other direction, a system can always be reduced to a single scalar equation of order $n$ by elimination, but the matrix route is the one that generalizes and the one that numerical software uses.

Example(A three-dimensional system with one real and one complex pair of eigenvalues)Original

Take the real system $x' = Ax$ with

A = \begin{pmatrix} 0 & 1 & 0 \\ -5 & -2 & 1 \\ 0 & 0 & -1 \end{pmatrix} , \qquad x_0 = \begin{pmatrix} 0 \\ 0 \\ 4\end{pmatrix} .

Written out, this is a damped oscillator $u'' + 2u' + 5u = w$ driven by a variable $w$ that decays on its own: $x = (u, u', w)^{\mathsf T}$ with $w' = -w$ .

1. Eigenvalues (step 1). The last row of $A$ has only the entry $-1$ on the diagonal, so expanding $\det(\lambda I_3 - A)$ along its last row gives

p(\lambda) = (\lambda + 1)\bigl(\lambda^2 + 2\lambda + 5\bigr) = (\lambda+1)\bigl((\lambda+1)^2 + 4\bigr) ,

using $\lambda^2 - (\operatorname{tr}B)\lambda + \det B$ with $B = \begin{pmatrix}0&1\\-5&-2\end{pmatrix}$ for the $2\times2$ block. The eigenvalues are $-1$ and $-1 \pm 2i$ : three distinct values, so $A$ is diagonalizable and step 5 is not needed. All three have real part $-1$ , so every solution decays like $e^{-t}$ , with oscillation of period $2\pi/2 = \pi$ .

2. The real eigenvalue (steps 2-3). Solve $(A + I_3)v = 0$ :

A + I_3 = \begin{pmatrix} 1 & 1 & 0 \\ -5 & -1 & 1 \\ 0 & 0 & 0\end{pmatrix} \quad\Longrightarrow\quad v_2 = -v_1, \quad v_3 = 5v_1 + v_2 = 4v_1 ,

the first equation giving $v_2 = -v_1$ and the second then giving $v_3$ . Taking $v_1 = 1$ ,

v = (1, -1, 4)^{\mathsf T}, \qquad \text{solution } \ e^{-t}(1,-1,4)^{\mathsf T} .

3. The complex pair (step 4). Take $\lambda = -1 + 2i$ and solve $(A - \lambda I_3)v = 0$ . The third row reads $(-1 - \lambda)v_3 = 0$ with $-1-\lambda = -2i \neq 0$ , so $v_3 = 0$ ; the first row then reads $-\lambda v_1 + v_2 = 0$ , so $v_2 = \lambda v_1$ . (The second row is then automatically satisfied, as it must be, since $\lambda$ is an eigenvalue.) With $v_1 = 1$ ,

v = (1,\, -1 + 2i,\, 0)^{\mathsf T} = p + iq, \qquad p = (1,-1,0)^{\mathsf T}, \quad q = (0,2,0)^{\mathsf T} ,

so $\alpha = -1$ , $\beta = 2$ and step 4 gives the two real solutions

e^{-t}\begin{pmatrix}\cos 2t \\ -\cos 2t - 2\sin 2t \\ 0\end{pmatrix} , \qquad e^{-t}\begin{pmatrix}\sin 2t \\ -\sin 2t + 2\cos 2t \\ 0\end{pmatrix} .

4. Fundamental matrix and the initial condition (steps 6-7). Collecting the three solutions as columns,

\Phi(t) = e^{-t}\begin{pmatrix} 1 & \cos 2t & \sin 2t \\ -1 & -\cos 2t - 2\sin 2t & 2\cos 2t - \sin 2t \\ 4 & 0 & 0 \end{pmatrix} , \qquad \Phi(0) = \begin{pmatrix} 1 & 1 & 0 \\ -1 & -1 & 2 \\ 4 & 0 & 0\end{pmatrix} .

Expanding $\det \Phi(0)$ along its last row gives $4\,(1\cdot 2 - 0\cdot(-1)) = 8 \neq 0$ , so $\Phi$ is a fundamental matrix. Solving $\Phi(0)c = (0,0,4)^{\mathsf T}$ : the third equation gives $4c_1 = 4$ , so $c_1 = 1$ ; the first gives $1 + c_2 = 0$ , so $c_2 = -1$ ; the second gives $-1 + 1 + 2c_3 = 0$ , so $c_3 = 0$ . Hence

x(t) = \Phi(t)\begin{pmatrix}1\\-1\\0\end{pmatrix} = e^{-t}\begin{pmatrix} 1 - \cos 2t \\ \cos 2t + 2\sin 2t - 1 \\ 4\end{pmatrix} .

5. Checks (step 8). At $t=0$ the bracket is $(1-1,\, 1 + 0 - 1,\, 4) = (0,0,4)$ , as required. The middle component must be the derivative of the first, since $x_2 = x_1'$ in this system:

\frac{d}{dt}\Bigl[e^{-t}(1 - \cos 2t)\Bigr] = e^{-t}\bigl(-(1-\cos 2t) + 2 \sin 2t\bigr) = e^{-t}(\cos 2t + 2\sin 2t - 1) ,

which it is. And the second equation of the system — row 2 of $A$ — namely $x_2' = -5x_1 - 2x_2 + x_3$ , holds: differentiating $x_2$ gives $e^{-t}(1 + 3\cos 2t - 4 \sin 2t)$ , while

-5e^{-t}(1-\cos 2t) - 2e^{-t}(\cos 2t + 2 \sin 2t - 1) + 4e^{-t} = e^{-t}\bigl((-5 + 2 + 4) + (5-2)\cos 2t - 4 \sin 2t\bigr) ,

the same thing. Finally $\operatorname{tr}A = 0 + (-2) + (-1) = -3$ , and indeed $\det\Phi(t) = e^{-3t}\det\Phi(0) = 8e^{-3t}$ , because each of the three columns carries one factor $e^{-t}$ and the remaining determinant is the constant $8$ (its entries involve $\cos 2t$ and $\sin 2t$ , but expanding along the last row leaves $4\bigl[\cos 2t(2\cos 2t - \sin 2t) - \sin 2t(-\cos 2t - 2 \sin 2t)\bigr] = 4\cdot 2 = 8$ by $\cos^2 + \sin^2 = 1$ ).

6. A number. At $t = \pi/4$ we have $\cos(\pi/2) = 0$ and $\sin(\pi/2) = 1$ , so

x(\pi/4) = e^{-\pi/4}\,(1,\ 1,\ 4)^{\mathsf T} \approx (0.4559381,\ 0.4559381,\ 1.8237525)^{\mathsf T} .

The displacement and the velocity happen to coincide a quarter of the way into the first oscillation.

7. The exponential, if it is wanted. $\exp(tA) = \Phi(t)\Phi(0)^{-1}$ , and $\Phi(0)^{-1} = \tfrac18\begin{pmatrix}0 & 0 & 2 \\ 8 & 0 & -2 \\ 4 & 4 & 0\end{pmatrix}$ , as one checks by multiplying out. Note that the solution above was obtained without ever forming this product: eigenvectors answer the initial-value problem directly, and the exponential is needed only when the propagator itself is the object of interest — as it is in the forced problem of nonhomogeneous systems.

ApplicationOriginal

Small oscillations of a mechanical or electrical network. A conservative system linearized about equilibrium obeys $M\ddot q + Kq = 0$ with $M$ positive definite and $K$ symmetric; written as a first-order system in $(q, \dot q)$ it is $x' = Ax$ with $A = \begin{pmatrix}0 & I_n \\ -M^{-1}K & 0\end{pmatrix}$ . Its eigenvalues satisfy $\lambda^2 = -\omega^2$ , where $\omega^2$ runs over the eigenvalues of $M^{-1}K$ (all real, since $M^{-1}K$ is self-adjoint for the inner product $\langle u, v\rangle_M = u^{\mathsf T}Mv$ ). At a stable equilibrium $K$ is positive definite as well as symmetric; then every $\omega^2 > 0$ , each $\pm i\omega$ is a conjugate pair on the imaginary axis, and step 4 of the method produces exactly the normal modes of normal modes, one mode per pair, oscillating at its own frequency $\omega$ . The two degenerate possibilities a merely symmetric $K$ allows are worth naming: $\omega^2 = 0$ gives a double eigenvalue $\lambda = 0$ with a defective block, hence a solution drifting linearly in $t$ (a zero mode, such as a rigid translation), and $\omega^2 < 0$ gives the real pair $\lambda = \pm|\omega|$ , an unstable direction growing exponentially rather than a mode. Adding linear damping moves an imaginary pair off the axis to the left, and its real part is the decay rate of that mode.
Stability of a steady state of a nonlinear system. Near an equilibrium $x_*$ of $x' = f(x)$ the deviation obeys $\xi' = Df(x_*)\xi$ to first order. The eigenvalues of the Jacobian therefore decide local stability, by the reading given at the end of the method, and this is the content of linearization: a sink when all real parts are negative, a saddle when they have both signs, oscillatory approach when there are complex pairs. The borderline cases — a purely imaginary pair, a zero eigenvalue — are precisely the ones the linear system cannot settle, and are where bifurcations occur.
Relaxation in coupled compartments. First-order chemical kinetics, radioactive decay chains, heat exchange between lumped bodies and linear pharmacokinetic models all have the form $c' = Ac$ with $A$ a matrix of rate constants whose columns sum to zero when the total amount is conserved. That column sum makes $(1,\dots,1)$ a left null vector, so $0$ is an eigenvalue and the corresponding eigenvector is the equilibrium distribution; the remaining eigenvalues, necessarily with negative real parts for a stable network, are the relaxation rates measured in a kinetics experiment. The same matrix structure, transposed, is the generator of a continuous-time Markov chain: with $Q = A^{\mathsf T}$ , which has zero row sums, $\exp(tQ) = \exp(tA)^{\mathsf T}$ is the matrix of transition probabilities over time $t$ .
Why software solves systems, not high-order equations. Every numerical integrator in use — Runge-Kutta, multistep, implicit — is written for first-order systems, so the reduction of the remark above is the first thing done to an $n$ -th order equation before it is handed to an ODE integrator. The eigenvalues of $A$ also set the step size an explicit method can take: the product of the step with the largest $|\lambda|$ must stay inside the method's stability region, which is why a system with widely separated eigenvalues is expensive to integrate explicitly.

Method

Nonhomogeneous systems and Duhamel's principle

Adapted

Add to PDF

Gives the solution of a forced linear system as the free evolution of the initial state plus a convolution of the forcing with the propagator, namely the Duhamel integral x(t) = exp(tA)x_0 plus the integral of exp((t-s)A)g(s) ds for constant A, and its fundamental-matrix form for time-dependent A, together with a recipe, the resonance that appears when the forcing matches an eigenvalue, and the reason the same formula reappears for the heat and wave equations.

Prerequisites: Linear systems and the matrix exponential, Functions of matrices and the matrix exponential, Variation of parameters, Determinants

Method(The structure of the solution set)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.5, p. 131

Let $A(t)$ and $g(t)$ be continuous on an interval $I$ and consider the nonhomogeneous (forced, inhomogeneous) system

x' = A(t)\,x + g(t) ,

whose homogeneous equation is the system $x' = A(t)x$ obtained by dropping $g$ .

If $z$ is any one solution of the nonhomogeneous system and $y$ any solution of the homogeneous one, then $y + z$ solves the nonhomogeneous system.
Conversely every solution of the nonhomogeneous system is of that form: if $z_1, z_2$ both solve it then $z_1 - z_2$ solves the homogeneous one.

So with $\Phi$ a fundamental matrix of the homogeneous system (linear systems) the general solution is

x(t) = \Phi(t)c + z(t), \qquad c \in \mathbb{C}^n ,

and the whole problem reduces to producing one particular solution $z$ . Hirsch, Smale and Devaney record both the reduction and the warning that goes with it: a particular solution can often be guessed, by undetermined coefficients, but guessing does not always work, whereas the variation-of-parameters formula below always does — subject only to one's being able to evaluate the integral it contains ((Hirsch, Smale & Devaney (2013), §6.5, p. 131)).

Idea of the proofOriginal

The map $x \mapsto x' - A(t)x$ is linear, and the two parts are the two halves of what linearity says about a linear equation: the difference of two solutions of $Lx = g$ solves $Lx = 0$ , and adding a solution of $Lx = 0$ to a solution of $Lx = g$ changes neither. Nothing about the particular form of $A$ or $g$ enters.

ProofOriginal

Write $L x = x' - A(t)x$ , defined on differentiable functions $I \to \mathbb{C}^n$ . $L$ is linear: $L(\alpha x + \beta y) = \alpha x' + \beta y' - A(\alpha x + \beta y) = \alpha Lx + \beta Ly$ , since differentiation is linear and $A(t)$ acts linearly at each $t$ .

If $Lz = g$ and $Ly = 0$ then $L(y+z) = Ly + Lz = g$ .
If $Lz_1 = Lz_2 = g$ then $L(z_1 - z_2) = g - g = 0$ .

For the displayed general solution, let $x$ solve the nonhomogeneous system. By part 2, $x - z$ solves the homogeneous one, so by part 3 of the existence statement in linear systems it equals $\Phi(t)c$ for exactly one $c$ ; conversely every such $\Phi(t)c + z(t)$ is a solution by part 1. $\square$

Method(Variation of constants: the Duhamel integral)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.5, Theorem (Variation of Parameters), p. 131

Let $A$ be a constant $n \times n$ matrix and let $g$ be continuous on an interval containing $0$ . Then

x(t) = \exp(tA)\Bigl[\, x_0 + \int_0^t \exp(-sA)\,g(s)\,ds \,\Bigr]

is a solution of $x' = Ax + g(t)$ with $x(0) = x_0$ ((Hirsch, Smale & Devaney (2013), §6.5, Theorem, p. 131)). Because $tA$ and $-sA$ commute, the bracket can be multiplied out to give the Duhamel integral form

x(t) = \underbrace{\exp(tA)\,x_0}_{\text{free evolution of the initial state}} \;+\; \underbrace{\int_0^t \exp\bigl((t-s)A\bigr)\, g(s)\, ds}_{\text{accumulated response to the forcing}} ,

and, this solution being the only one with $x(0) = x_0$ (part 1 of the existence statement in linear systems applies to the nonhomogeneous system as well), the two displays are the same function.

Read the second display pointwise in $s$ : the forcing $g(s)\,ds$ delivered during $[s, s+ds]$ acts as a fresh initial condition at time $s$ , and is then propagated to time $t$ by the propagator $\exp((t-s)A)$ of the unforced system. The total response is the superposition of these, and it is a convolution of the forcing with the propagator: only the elapsed time $t - s$ appears, and only values $s \le t$ contribute, so the response is causal. The matrix

K(t,s) = \exp\bigl((t-s)A\bigr) \quad (s \le t), \qquad K(t,s) = 0 \quad (s > t)

is the impulse response, or causal Green's function, of the system (causal Green's functions).

Idea of the proofOriginal

Freeze the forcing and the answer is $\exp(tA)c$ with $c$ constant. Switch the forcing on, and the only way the formula can fail is for $c$ to have to drift: so look for a solution $x(t) = \exp(tA)\,u(t)$ with $u$ an unknown function — this is what "varying the constant" means. Substituting, the two terms produced by differentiating the product are $A\exp(tA)u$ , which is exactly $Ax$ and cancels against the left-hand side, and $\exp(tA)u'$ , which must therefore equal $g$ . That leaves $u' = \exp(-tA)g$ , a formula to be integrated and nothing more. The structure survives verbatim for time-dependent $A$ , with a fundamental matrix in place of the exponential, because the cancellation used only $\Phi' = A\Phi$ .

ProofAdaptedAdapted from Hirsch, Smale & Devaney (2013), §6.5, p. 132

The integral is differentiable. $s \mapsto \exp(-sA)g(s)$ is continuous, being a product of a continuous matrix-valued function and a continuous vector-valued one. So by the fundamental theorem of calculus in each of the $n$ components separately, $t \mapsto \int_0^t \exp(-sA)g(s)\,ds$ is differentiable with derivative $\exp(-tA)g(t)$ . Continuity of $g$ is used here and nowhere else.
Differentiate the product. Write $v(t) = x_0 + \int_0^t \exp(-sA)g(s)\,ds$ , so that $x = \exp(tA)v$ . The product rule for a matrix-valued and a vector-valued factor, together with $\frac{d}{dt}\exp(tA) = A\exp(tA)$ (matrix exponential) and step 1, gives

x'(t) = A \exp(tA)\,v(t) + \exp(tA)\exp(-tA)\,g(t) = A\,x(t) + g(t) ,

the second term collapsing because $\exp(tA)\exp(-tA) = I$ . This is the route of (Hirsch, Smale & Devaney (2013), §6.5, p. 132), which differentiates the formula as it stands rather than deriving it.

The initial value. At $t = 0$ the integral is empty and $\exp(0) = I$ , so $x(0) = x_0$ .
The two displays agree. For each fixed $s$ the matrices $tA$ and $-sA$ commute, so the product rule for exponentials of commuting matrices gives $\exp(tA)\exp(-sA) = \exp((t-s)A)$ . Since $\exp(tA)$ is a constant matrix as far as the $s$ -integration is concerned, it may be taken inside the integral:

\exp(tA)\int_0^t \exp(-sA)g(s)\,ds = \int_0^t \exp(tA)\exp(-sA)g(s)\,ds = \int_0^t \exp\bigl((t-s)A\bigr)g(s)\,ds .

Uniqueness. If $x$ and $\tilde x$ both solve the problem, their difference solves the homogeneous system with initial value $0$ , hence vanishes identically (linear systems). So the formula gives the solution, not merely one. $\square$

Method(Time-dependent coefficients)Original

Let $A(t)$ and $g(t)$ be continuous on an interval $I$ containing $t_0$ , and let $\Phi$ be any fundamental matrix of $x' = A(t)x$ on $I$ . Then the unique solution of

x' = A(t)x + g(t), \qquad x(t_0) = x_0

on $I$ is

x(t) = \Phi(t)\Phi(t_0)^{-1}x_0 + \int_{t_0}^{t} \Phi(t)\,\Phi(s)^{-1}\,g(s)\,ds ,

and the answer does not depend on which fundamental matrix was used. Writing $U(t,s) = \Phi(t)\Phi(s)^{-1}$ for the propagator, this is again

x(t) = U(t,t_0)\,x_0 + \int_{t_0}^t U(t,s)\,g(s)\,ds ,

identical in shape to the constant-coefficient formula, of which it is the case $U(t,s) = \exp((t-s)A)$ . What is lost for non-constant $A$ is only that $U(t,s)$ need no longer depend on $t - s$ alone, and that there is in general no closed form for it: unlike the constant-coefficient case, $\Phi(t) = \exp\bigl(\int_{t_0}^t A\bigr)$ is false in general, because $A(t)$ and $\int A$ need not commute. Neither cited source states this version; it is ours.

Idea of the proofOriginal

Run the previous idea with $\Phi$ in place of $\exp(tA)$ . Substituting $x = \Phi u$ and using $\Phi' = A\Phi$ makes the $Ax$ terms cancel and leaves $\Phi u' = g$ , which can be solved for $u'$ precisely because a fundamental matrix is invertible at every point — the one fact about fundamental matrices that does all the work. Integrating $u'$ from $t_0$ and fitting $u(t_0)$ to the initial condition gives the formula. Independence of the choice of $\Phi$ then needs no computation: two fundamental matrices differ by a constant nonsingular factor on the right, which cancels in every $\Phi(t)\Phi(s)^{-1}$ .

ProofOriginal

$\Phi^{-1}$ is continuous. By part 2 of the existence statement in linear systems, $\det \Phi(s) \neq 0$ for every $s \in I$ . The entries of $\Phi(s)^{-1}$ are cofactors of $\Phi(s)$ divided by $\det\Phi(s)$ (determinants), hence continuous in $s$ , being polynomials in continuous functions divided by a nonvanishing continuous function. So $s \mapsto \Phi(s)^{-1}g(s)$ is continuous and its integral is differentiable, with derivative $\Phi(t)^{-1}g(t)$ , by the fundamental theorem of calculus.
The substitution. Put

u(t) = \Phi(t_0)^{-1}x_0 + \int_{t_0}^{t}\Phi(s)^{-1}g(s)\,ds, \qquad x(t) = \Phi(t)u(t) ,

which is the displayed formula once $\Phi(t)$ is taken inside the integral, legitimate because $\Phi(t)$ is constant in $s$ . Differentiating the product and using step 1 and $\Phi' = A\Phi$ ,

x' = \Phi'(t)u(t) + \Phi(t)\Phi(t)^{-1}g(t) = A(t)\Phi(t)u(t) + g(t) = A(t)x(t) + g(t) .

The initial value. The integral vanishes at $t = t_0$ , so $x(t_0) = \Phi(t_0)\Phi(t_0)^{-1}x_0 = x_0$ .
Uniqueness. Two solutions differ by a solution of the homogeneous system vanishing at $t_0$ , which is identically zero.
Independence of $\Phi$ . If $\Psi$ is another fundamental matrix then $\Psi = \Phi Q$ for a constant nonsingular $Q$ (part 3 of the existence statement in linear systems), so for all $t, s$

\Psi(t)\Psi(s)^{-1} = \Phi(t)QQ^{-1}\Phi(s)^{-1} = \Phi(t)\Phi(s)^{-1} .

Both terms of the formula are built only from such products, so neither changes.

Why $\exp\bigl(\int A\bigr)$ fails. Set $B(t) = \int_{t_0}^t A(\tau)\,d\tau$ . Differentiating the series for $\exp(B(t))$ term by term would need $\frac{d}{dt}B(t)^k = \sum_{j} B^{j}B'B^{k-1-j}$ to collapse to $kB^{k-1}B'$ , which requires $B(t)$ and $B'(t) = A(t)$ to commute. They need not: for

A(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix} + t\begin{pmatrix}0&0\\1&0\end{pmatrix}

the two constituent matrices do not commute, and $A(t)$ does not commute with $\int_0^t A$ . (This is the same obstruction that makes $\exp(X)\exp(Y) \ne \exp(X+Y)$ for non-commuting $X, Y$ ; see the matrix exponential.) $\square$

Remark(Duhamel's principle)AdaptedAdapted from Evans (2010), §2.3.1(c), p. 47; Strauss (2008), §3.4, pp. 77-78

The formula above is one instance of a principle with a name. Evans states it for the heat equation — the solution of $u_t - \Delta u = f$ with zero initial data is built by integrating, over the switch-on time $s$ , the solutions of the homogeneous problem started at time $s$ with initial data $f(\cdot, s)$ — and records in a footnote that the principle applies widely to linear ordinary and partial differential equations and does not depend on the particular structure of the heat equation ((Evans (2010), §2.3.1(c), p. 47)). Strauss puts the same thing operationally: a homogeneous problem one already knows how to solve delivers the forced problem as well, and that implication is what the name Duhamel's principle covers ((Strauss (2008), §3.4, p. 78)).

Strauss's own ODE template is second order rather than first:

\frac{d^2u}{dt^2} + A^2 u = f(t), \quad u(0) = \varphi, \quad u'(0) = \psi \qquad\Longrightarrow\qquad u(t) = S'(t)\varphi + S(t)\psi + \int_0^t S(t-s)f(s)\,ds ,

where $S(t) = A^{-1}\sin tA$ and $S'(t) = \cos tA$ , and $A^2$ may be a positive constant or a positive square matrix ((Strauss (2008), §3.4, (13)-(15), p. 77)). The key, as he puts it, is that $S(t)\psi$ is the solution of the problem with $\varphi = 0$ and $f = 0$ : the solution operator attached to the velocity datum is the one that also convolves against the forcing. Our first-order formula is the same statement with $U(t,s)$ in that role, and the two are consistent — writing the second-order equation as a first-order system in $(u, u')$ turns $U(t,0)$ into the block matrix $\begin{pmatrix} \cos tA & A^{-1}\sin tA \\ -A \sin tA & \cos tA\end{pmatrix}$ , whose top-right block is $S(t)$ .

The same pattern, with the same name, produces the solution of the forced wave equation and of the forced heat equation from their homogeneous solution operators; those are taken up in Green's functions for the wave and heat equations. The abstract version, in which $\exp(tA)$ is replaced by a strongly continuous semigroup on a Banach space, is in operator semigroups.

Remark(Comparison with scalar variation of parameters, and what the hypotheses do)AdaptedAdapted from Coddington (1989), Ch. 2, §10, Theorem 20, (10.3) and (10.4), pp. 85-86

The scalar case is the same calculation. For $L(y) = y^{(n)} + a_1 y^{(n-1)} + \dots + a_n y = b$ , Coddington writes a particular solution as $\psi_p = \sum_k u_k \phi_k$ with $\phi_1, \dots, \phi_n$ a basis of solutions of $L(y) = 0$ , and imposes the $n$ conditions

\sum_k u_k' \phi_k^{(j)} = 0 \ \ (j = 0, \dots, n-2), \qquad \sum_k u_k'\phi_k^{(n-1)} = b ,

a linear system whose coefficient determinant is the Wronskian $W(\phi_1, \dots, \phi_n)$ and is therefore nonzero; solving it by Cramer's rule and integrating from $x_0$ gives his formula (10.3), and he records that the resulting $\psi_p$ satisfies $\psi_p(x_0) = \psi_p'(x_0) = \dots = \psi_p^{(n-1)}(x_0) = 0$ ((Coddington (1989), Ch. 2, §10, Theorem 20, (10.3)-(10.4), pp. 85-86)). Under the reduction of a scalar equation to a system (linear systems), that system of conditions is exactly $\Phi u' = g$ with $g = (0,\dots,0,b)^{\mathsf T}$ , the coefficient determinant is $\det \Phi$ , and the vanishing initial data are the statement $x(t_0) = 0$ for the Duhamel term. The matrix derivation is the same argument with the bookkeeping done once instead of $n$ times; see also variation of parameters.

Where the hypotheses are used.

Continuity of $g$ is used only in step 1 of each proof, to differentiate the integral. If $g$ is merely piecewise continuous — a switched-on forcing, say — the formula still defines a continuous function, it still satisfies the differential equation at every continuity point of $g$ , and $x$ is then continuously differentiable except at the jumps of $g$ , where $x'$ jumps by the jump of $g$ . Coddington makes the same observation for the scalar formula ((Coddington (1989), Ch. 2, §10, Exercise 5, p. 90)). An impulsive forcing is outside the statement altogether and is handled by a jump in the initial condition, or by the transform method of Laplace transforms for initial-value problems.
$A$ constant is what makes $U(t,s)$ depend on $t-s$ , hence makes the Duhamel term a convolution; for variable $A$ it is still an integral operator, but not a convolution.
Invertibility of $\Phi$ at every point is what the substitution $x = \Phi u$ needs, and it is not an extra assumption: it is automatic for a fundamental matrix.
No hypothesis whatever is needed on the size of $g$ : the formula is exact, not an approximation, and there is no smallness or resonance condition. What resonance does is make the integral large, not invalid — see part 2 of the example.

Method(Solving x' = Ax + g(t))Original

Solve the homogeneous system first. Get $\exp(tA)$ , or a fundamental matrix $\Phi$ , by the recipe in linear systems. Nothing about the forcing is needed yet.
Decompose the forcing in the eigenbasis, when there is one. With $A$ diagonalizable, write $g(t) = \sum_j \gamma_j(t)\,v_j$ with $Av_j = \lambda_j v_j$ ; then $\exp((t-s)A)g(s) = \sum_j e^{\lambda_j(t-s)}\gamma_j(s)\,v_j$ and the $n$ -dimensional integral splits into $n$ scalar ones,

x(t) = \sum_j \Bigl[ c_j e^{\lambda_j t} + \int_0^t e^{\lambda_j (t-s)}\gamma_j(s)\,ds \Bigr] v_j .

This is almost always the cheapest route by hand, and it shows at once which modes the forcing excites: a $\gamma_j$ that vanishes leaves mode $j$ unforced. 3. Or integrate the Duhamel formula directly, $x(t) = \exp(tA)x_0 + \int_0^t \exp((t-s)A)g(s)\,ds$ , when $\exp(tA)$ is available in closed form. 4. Or guess, when the forcing is of exponential-polynomial type. If $g(t) = e^{\mu t}w$ with $\mu$ not an eigenvalue of $A$ , try $z(t) = e^{\mu t}\zeta$ : substituting gives $(\mu I - A)\zeta = w$ , solvable because $\mu I - A$ is then invertible, so $\zeta = (\mu I - A)^{-1}w$ . For $g(t) = e^{\mu t}\bigl(w_0 + tw_1 + \dots\bigr)$ take a polynomial ansatz of the same degree, raising the degree by the size of the largest Jordan block of $\mu$ if $\mu$ is an eigenvalue. A sinusoidal real forcing $g(t) = \operatorname{Re}(e^{i\omega t}w)$ is the case $\mu = i\omega$ , and the transfer matrix $(i\omega I - A)^{-1}$ is the object a physicist calls the system's frequency response. 5. Add the free term and fit the initial condition. $x(t) = \exp(tA)x_0 + z(t) - \exp(tA)z(0)$ if $z$ is a particular solution obtained by guessing; or read $x_0$ straight off the Duhamel form, where the initial condition is already built in. 6. Checks. $x(0) = x_0$ ; $x'(0) = Ax_0 + g(0)$ ; and if $g$ is constant and $A$ invertible, $x(t) \to -A^{-1}g$ as $t \to \infty$ whenever every eigenvalue has negative real part, since then the steady state solves $0 = Ax_\infty + g$ .

What to expect. If every eigenvalue has $\operatorname{Re}\lambda < 0$ , the term $\exp(tA)x_0$ is a transient that dies out and the Duhamel term is the forced response that survives; the long-time behaviour is then a property of the forcing and not of the initial condition. If the forcing is periodic there is exactly one periodic solution and everything converges to it (see the application). If $\mu$ coincides with an eigenvalue, the ansatz of step 4 fails and the response grows by an extra factor $t$ per repetition: that is resonance.

Example(Two coupled compartments with a source)Original

Two well-stirred compartments exchange material with each other at rate $1$ and each loses material to the outside at rate $2$ ; a source feeds compartment $1$ . With $c = (c_1, c_2)^{\mathsf T}$ the amounts present,

c_1' = -2c_1 - (c_1 - c_2) + f(t), \qquad c_2' = -2c_2 - (c_2 - c_1) ,

that is $c' = Ac + g(t)$ with

A = \begin{pmatrix} -3 & 1 \\ 1 & -3\end{pmatrix}, \qquad g(t) = \begin{pmatrix} f(t) \\ 0\end{pmatrix} , \qquad c(0) = \begin{pmatrix}0\\0\end{pmatrix} .

1. The homogeneous system (step 1). $A = -3I + \begin{pmatrix}0&1\\1&0\end{pmatrix}$ , and the swap matrix has eigenvectors $v_+ = (1,1)^{\mathsf T}$ and $v_- = (1,-1)^{\mathsf T}$ with eigenvalues $+1$ and $-1$ . Adding $-3$ to each,

Av_+ = -2v_+, \qquad Av_- = -4 v_- :

$\lambda_+ = -2$ is the slow mode (both compartments together, decaying only through the loss to the outside) and $\lambda_- = -4$ the fast one (the imbalance between them, which also relaxes through the exchange). Both are negative, so the unforced system decays.

2. Decompose the forcing (step 2). Since $(1,0)^{\mathsf T} = \tfrac12 v_+ + \tfrac12 v_-$ , the forcing is $g(t) = \tfrac12 f(t) v_+ + \tfrac12 f(t) v_-$ : it drives both modes equally.

Part 1: a constant source, $f \equiv 1$ . Step 2's formula with $c(0) = 0$ gives

c(t) = \tfrac12\Bigl[\int_0^t e^{-2(t-s)}ds\Bigr]v_+ + \tfrac12\Bigl[\int_0^t e^{-4(t-s)}ds\Bigr]v_- = \frac{1 - e^{-2t}}{4}\,v_+ + \frac{1-e^{-4t}}{8}\,v_- ,

using $\int_0^t e^{-a(t-s)}ds = (1 - e^{-at})/a$ . Componentwise,

c_1(t) = \frac38 - \frac14 e^{-2t} - \frac18 e^{-4t} , \qquad c_2(t) = \frac18 - \frac14 e^{-2t} + \frac18 e^{-4t} .

Checks. At $t = 0$ : $\tfrac38 - \tfrac14 - \tfrac18 = 0$ and $\tfrac18 - \tfrac14 + \tfrac18 = 0$ . As $t \to \infty$ the limit is $(\tfrac38, \tfrac18)^{\mathsf T}$ , and indeed $-A^{-1}g = \tfrac18\begin{pmatrix}3&1\\1&3\end{pmatrix}\begin{pmatrix}1\\0\end{pmatrix} = (\tfrac38, \tfrac18)^{\mathsf T}$ , since $\det A = 8$ and $A^{-1} = \tfrac18\begin{pmatrix}-3&-1\\-1&-3\end{pmatrix}$ . Substituting into the first equation, $c_1' = \tfrac12 e^{-2t} + \tfrac12 e^{-4t}$ while

-3c_1 + c_2 + 1 = \Bigl(-\tfrac98 + \tfrac18 + 1\Bigr) + \Bigl(\tfrac34 - \tfrac14\Bigr)e^{-2t} + \Bigl(\tfrac38 + \tfrac18\Bigr)e^{-4t} = \tfrac12 e^{-2t} + \tfrac12 e^{-4t} ,

the same. A number. At $t = \tfrac12\ln 2$ , where $e^{-2t} = \tfrac12$ and $e^{-4t} = \tfrac14$ ,

c_1 = \tfrac38 - \tfrac18 - \tfrac1{32} = \tfrac{7}{32} = 0.21875 , \qquad c_2 = \tfrac18 - \tfrac18 + \tfrac1{32} = \tfrac1{32} = 0.03125 ,

so the second compartment holds exactly one seventh of the first at that moment.

Part 2: a source that resonates with the slow mode, $f(t) = e^{-2t}$ . Now $\mu = -2$ is an eigenvalue, so the guess of step 4 is unavailable. The Duhamel integrals of step 2 are still elementary: the $v_+$ integral has $e^{-2(t-s)}e^{-2s} = e^{-2t}$ , with no $s$ left in it, so it contributes $t e^{-2t}$ , while the $v_-$ integral is an ordinary exponential one:

c(t) = \tfrac12 \Bigl[\int_0^t e^{-2t}\,ds\Bigr]v_+ + \tfrac12\Bigl[\int_0^t e^{-4(t-s)}e^{-2s}ds\Bigr]v_- = \frac{t}{2}e^{-2t}\,v_+ + \frac{e^{-2t} - e^{-4t}}{4}\,v_- ,

the second integral being $e^{-4t}\int_0^t e^{2s}ds = e^{-4t}(e^{2t}-1)/2$ . Componentwise,

c_1(t) = \frac t2 e^{-2t} + \frac{e^{-2t}-e^{-4t}}{4} , \qquad c_2(t) = \frac t2 e^{-2t} - \frac{e^{-2t}-e^{-4t}}{4} .

Checks. Both vanish at $t = 0$ . Differentiating $c_1$ gives $\bigl(\tfrac12 - t\bigr)e^{-2t} - \tfrac12 e^{-2t} + e^{-4t} = -te^{-2t} + e^{-4t}$ , while

-3c_1 + c_2 + e^{-2t} = -3\Bigl(\tfrac t2 e^{-2t}\Bigr) + \tfrac t2 e^{-2t} - \tfrac{3+1}{4}\bigl(e^{-2t}-e^{-4t}\bigr) + e^{-2t} = -te^{-2t} + e^{-4t} ,

the same. A number. At $t = \tfrac12\ln 2$ ,

c_1 = \frac{\ln 2}{8} + \frac1{16} = 0.149143398 , \qquad c_2 = \frac{\ln 2}{8} - \frac1{16} = 0.024143398 .

3. What resonance did and did not do. The secular factor $t$ appears only in the $v_+$ component, the mode whose eigenvalue the forcing matched; the $v_-$ component is an ordinary combination of exponentials. And because $\operatorname{Re}\lambda_+ = -2 < 0$ , the product $te^{-2t}$ still decays: resonance with a decaying mode produces a bounded bump, peaking at $t = 1/2$ in the $v_+$ coefficient, not unbounded growth. Unbounded growth needs the resonant eigenvalue to sit on the imaginary axis, as it does for an undamped oscillator driven at its natural frequency.

ApplicationOriginal

Transient and steady state in a driven circuit or oscillator. For a damped linear system driven at a single frequency, $g(t) = \operatorname{Re}(e^{i\omega t}w)$ , step 4 of the method gives the steady state $\operatorname{Re}\bigl(e^{i\omega t}(i\omega I - A)^{-1}w\bigr)$ , and the Duhamel formula shows that everything else — the dependence on how the system was switched on — is confined to $\exp(tA)x_0$ , which decays. The matrix $(i\omega I - A)^{-1}$ is the frequency response: its magnitude is the gain, its argument the phase lag, and its poles are at $i\omega = \lambda$ , which is why the response peaks near the imaginary parts of the eigenvalues. Hirsch, Smale and Devaney prove the corresponding exact statement for the damped oscillator forced by $\cos t$ : there is exactly one $2\pi$ -periodic solution, found by solving $\bigl(\exp(2\pi A) - I\bigr)x_0 = -W$ — solvable because every eigenvalue of $\exp(2\pi A)$ has modulus below $1$ when every eigenvalue of $A$ has negative real part, so $1$ is not among them — and every other solution converges to it ((Hirsch, Smale & Devaney (2013), §6.5, pp. 132-133)).
Impulse response and the measurement of a linear system. Because the response to $g$ is $\int_0^t K(t,s)g(s)\,ds$ with $K(t,s) = \exp((t-s)A)$ , the system is characterized completely by what it does to a sharp kick: feed in $g(s) \approx \delta(s)e_j$ and the output is the $j$ -th column of $\exp(tA)$ . This is why experimentalists measure impulse or step responses rather than solving for $A$ — and why, conversely, deconvolving a measured response against a known forcing recovers the kernel. The same kernel, extended to $s > t$ by zero, is the causal Green's function of causal Green's functions.
Heat conduction with a time-dependent source. Discretizing $u_t = \kappa u_{xx} + f(x,t)$ in space (finite differences) gives exactly $x' = Ax + g(t)$ with $A$ the discrete Laplacian, and the Duhamel formula is the discrete version of Evans's heat-equation formula: each slice of heat released at time $s$ is diffused for the remaining time $t - s$ . Running the same formula with $A$ replaced by the continuum Laplacian and $\exp((t-s)A)$ by the heat kernel gives the continuum statement (the heat kernel).
Exponential integrators in numerical practice. Over one step of length $h$ , the exact solution is $x(t+h) = \exp(hA)x(t) + \int_0^h \exp((h-\tau)A)g(t+\tau)\,d\tau$ . Approximating only the integral — for example by freezing $g$ at $g(t)$ , which gives $x(t+h) \approx \exp(hA)x(t) + A^{-1}(\exp(hA) - I)g(t)$ — yields a method that is exact for the stiff linear part and therefore stable at step sizes an explicit method could not use. This is the basis of exponential time differencing, and it is the practical answer to the stiffness described in stiffness.
Linear control and the variation-of-constants formula. With $g(t) = Bu(t)$ for a control input $u$ , the Duhamel formula says the reachable states at time $T$ from the origin are the values of $\int_0^T \exp((T-s)A)Bu(s)\,ds$ . The span of those values, as $u$ ranges over all inputs, is the column space of $\bigl(B, AB, \dots, A^{n-1}B\bigr)$ — the controllability criterion — and the formula above is where that criterion comes from.

Theorem · advanced

Floquet theory

Adapted

Add to PDF

For a linear system whose coefficient matrix has period T, the fundamental matrix normalised at the origin factors as a T-periodic matrix times the exponential of a constant one, so the system is a periodic change of variable away from a constant-coefficient system; the eigenvalues of the monodromy matrix — the Floquet multipliers — do not depend on any choice made and decide decay, boundedness and the existence of periodic solutions, which for Hill's equation reduces to the single test |tr M| < 2.

Prerequisites: Linear systems and the matrix exponential, Structure of linear equations and the Wronskian, Existence and uniqueness of solutions, Functions of matrices and the matrix exponential, Jordan normal form, Eigenvalues and eigenvectors, Determinants

Definition(Periodic systems, the monodromy matrix and the Floquet multipliers)AdaptedAdapted from Coddington (1989), Ch. 1, Exercise 6, p. 46, and Ch. 3, §3, Exercise 6, p. 110

Let $A$ be an $n \times n$ matrix of continuous complex-valued functions on all of $\mathbb{R}$ , and let $T > 0$ satisfy

A(t + T) = A(t) \qquad \text{for every } t \in \mathbb{R} .

We call $x' = A(t)x$ a $T$ -periodic linear system. Here $T$ is a period, not necessarily the smallest one; every statement below holds for whichever period is chosen, and the conclusions change with it, which is a point we return to.

Because $A$ is continuous on $\mathbb{R}$ , the initial-value problem for the system has a unique solution defined on the whole line (Linear systems and the matrix exponential). Write $\Phi$ for the principal fundamental matrix at $0$ : the matrix solution of

\Phi'(t) = A(t)\,\Phi(t), \qquad \Phi(0) = I_n ,

whose $j$ -th column is the solution with initial vector $e_j$ . The monodromy matrix of the system relative to this normalisation is

M = \Phi(T) ,

the matrix that advances the state by one period: the solution with $x(0) = x_0$ has $x(T) = Mx_0$ . Its eigenvalues $\rho_1, \dots, \rho_n$ , listed with algebraic multiplicity, are the Floquet multipliers (also called characteristic multipliers, (Strogatz (1994), §8.7, p. 282)) of the system. A Floquet exponent attached to $\rho_j$ is any complex number $\mu_j$ with

e^{\mu_j T} = \rho_j ,

so $\mu_j$ is determined only modulo $2\pi i/T$ , while $\operatorname{Re}\mu_j = T^{-1}\log|\rho_j|$ is determined exactly.

In dimension one the whole construction collapses to a single number: for $y' + a(t)y = 0$ with $a$ of period $T$ , the solution normalised by $y(0) = 1$ is $\exp\bigl(-\int_0^t a\bigr)$ , so $M$ is the $1 \times 1$ matrix $\exp\bigl(-\int_0^T a\bigr)$ and this is the only multiplier (Coddington (1989), Ch. 1, Exercise 6(b), p. 46). In dimension two, $M$ is the matrix of the four constants that express the translated basis $t \mapsto \varphi_k(t + T)$ in terms of the basis $\varphi_1, \varphi_2$ normalised at $0$ (Coddington (1989), Ch. 3, §3, Exercise 6(c)-(d), p. 110). Coddington poses both of these as exercises; the general statement below is not in the library, and the blocks that prove it are labelled accordingly.

Theorem(Floquet's theorem)Original

Let $x' = A(t)x$ be a $T$ -periodic linear system as above, with principal fundamental matrix $\Phi$ at $0$ and monodromy matrix $M = \Phi(T)$ .

Translation. $\Phi(t + T) = \Phi(t)\,M$ for every $t \in \mathbb{R}$ , and more generally $\Phi(t + kT) = \Phi(t)M^k$ for every $k \in \mathbb{Z}$ .
Nonsingularity. $M$ is invertible; indeed $\det M = \exp \int_0^{T} \operatorname{tr} A(s)\, ds ,$ so no multiplier is zero.
The multipliers are intrinsic. If $\Psi$ is any fundamental matrix of the system, then $\Psi(t + T) = \Psi(t)\,M_\Psi$ with $M_\Psi = \Psi(0)^{-1} M \Psi(0)$ . Monodromy matrices built from different fundamental matrices are therefore similar, and the multipliers $\rho_1, \dots, \rho_n$ depend only on $A$ and $T$ .
Normal form. There are a matrix $P$ of continuously differentiable functions on $\mathbb{R}$ , with $P(t+T) = P(t)$ and $\det P(t) \neq 0$ for every $t$ , and a constant matrix $R$ , such that $\Phi(t) = P(t)\, e^{tR}, \qquad e^{TR} = M, \qquad P(0) = I_n .$ The substitution $x = P(t)y$ turns $x' = A(t)x$ into the constant-coefficient system $y' = Ry$ . The eigenvalues of $R$ are a set of Floquet exponents, one for each multiplier.
Floquet solutions. For each multiplier $\rho$ there is a non-trivial solution with $x(t + T) = \rho\, x(t) \qquad \text{for all } t,$ equivalently $x(t) = e^{\mu t} p(t)$ with $p$ of period $T$ and $e^{\mu T} = \rho$ . Such a solution is called a Floquet solution, or a Bloch solution.

Theorem(What the multipliers decide)AdaptedAdapted from Coddington (1989), Ch. 3, §4, Exercise 3, p. 116; Strogatz (1994), §8.7, pp. 281-282

With the notation above, let $\rho_1, \dots, \rho_n$ be the Floquet multipliers and $\mu_1, \dots, \mu_n$ corresponding exponents.

Decay. Every solution satisfies $x(t) \to 0$ as $t \to +\infty$ exactly when every multiplier obeys $|\rho_j| < 1$ ; equivalently, when every exponent has $\operatorname{Re}\mu_j < 0$ . In that case the decay is exponential: for each $r$ with $\max_j |\rho_j| < r < 1$ there is a constant $C$ with $\|x(t)\| \le C\, r^{\,t/T}\|x(0)\|$ for all $t \ge 0$ . (This is the criterion Strogatz (Strogatz (1994), §8.7, p. 281) states for the non-trivial multipliers of a closed orbit of an autonomous system; the application below connects that setting to this one.)
Boundedness. Every solution is bounded on $[0, \infty)$ if and only if $|\rho_j| \le 1$ for every $j$ and every multiplier of modulus $1$ is a semisimple eigenvalue of $M$ , meaning that its algebraic and geometric multiplicities agree.
Instability. If $|\rho_j| > 1$ for some $j$ , there is a solution with $\|x(t)\| \to \infty$ as $t \to +\infty$ .
Periodic and antiperiodic solutions. A non-trivial solution of period $T$ exists precisely when $1$ occurs among the multipliers, and a non-trivial solution obeying $x(t + T) = -x(t)$ , hence of period $2T$ , precisely when $-1$ does.
Hill's equation. Let $\alpha$ be real, continuous and of period $T$ , and consider $y'' + \alpha(t)y = 0$ , written as the system $x = (y, y')$ with $A(t) = \begin{pmatrix} 0 & 1 \\ -\alpha(t) & 0\end{pmatrix}$ . Then $\operatorname{tr}A \equiv 0$ , so $\det M = 1$ and the multipliers are the roots of $\rho^2 - \tau\rho + 1 = 0$ with $\tau = \operatorname{tr} M$ . Consequently:
- if $|\tau| < 2$ , every solution is bounded on $\mathbb{R}$ ;
- if $|\tau| > 2$ , some solution grows exponentially as $t \to +\infty$ ;
- $\tau = 2$ exactly when there is a non-trivial solution of period $T$ , and $\tau = -2$ exactly when there is a non-trivial solution with $y(t+T) = -y(t)$ (Coddington (1989), Ch. 3, §4, Exercise 3(b)-(c), p. 116); in these two borderline cases all solutions are bounded if and only if $M = I_2$ , respectively $M = -I_2$ .
In the notation of Coddington's exercise (Coddington (1989), Ch. 3, §4, Exercise 3, p. 116), where $\varphi_1, \varphi_2$ are the solutions with $\varphi_1(0) = 1$ , $\varphi_1'(0) = 0$ and $\varphi_2(0) = 0$ , $\varphi_2'(0) = 1$ , the Wronskian $W(\varphi_1, \varphi_2)$ is identically $1$ and $\tau = \varphi_1(T) + \varphi_2'(T)$ .

Idea of the proofOriginal

Everything follows from one observation: shifting time by a period carries solutions to solutions, because the coefficient matrix is unchanged by the shift. The shift is therefore a linear map of the $n$ -dimensional solution space to itself, and $M$ is its matrix in the basis normalised at $t = 0$ . So a global question about a system with variable coefficients — do solutions grow? — becomes a finite-dimensional question about the powers of one matrix.

The normal form says more: since $M$ is invertible it is the exponential of some matrix $B$ (build $B$ from the Jordan form of $M$ , using a logarithm of each eigenvalue and the terminating logarithmic series on the nilpotent part), and dividing $B$ by $T$ gives a constant $R$ whose exponential advances exactly like $M$ over one period. Stripping $e^{tR}$ off $\Phi$ leaves something that repeats with period $T$ . The growth of solutions is carried entirely by $e^{tR}$ , and the periodic factor only wobbles.

ProofOriginal

Step 1: translation. Put $\Theta(t) = \Phi(t + T)$ . Then $\Theta'(t) = \Phi'(t+T) = A(t+T)\Phi(t+T) = A(t)\Theta(t)$ , where the last equality is the periodicity of $A$ — this is the only place it is used. So $\Theta$ solves the same matrix equation as $\Phi$ and has $\Theta(0) = \Phi(T) = M$ . The matrix $t \mapsto \Phi(t)M$ solves it too, since $(\Phi M)' = \Phi'M = A\Phi M$ , and takes the same value $M$ at $t = 0$ . By uniqueness for the initial-value problem (Linear systems and the matrix exponential) applied to each column, $\Phi(t+T) = \Phi(t)M$ for all $t$ . Iterating gives $\Phi(t + kT) = \Phi(t)M^k$ for $k \ge 0$ , and the relation $\Phi(t) = \Phi(t - T)M$ gives it for negative $k$ once $M$ is known to be invertible, which is Step 2.

Step 2: the determinant. Let $w(t) = \det \Phi(t)$ . The determinant is linear in each row separately, so differentiating it row by row,

w'(t) = \sum_{j=1}^{n} \det \Phi_{[j]}(t) ,

where $\Phi_{[j]}$ is $\Phi$ with its $j$ -th row replaced by the derivative of that row. Row $j$ of $\Phi' = A\Phi$ is $\sum_{k} a_{jk}(t)\,(\text{row } k \text{ of } \Phi)$ . Expanding $\det\Phi_{[j]}$ by linearity in that row, every term with $k \neq j$ is a determinant with two equal rows and vanishes, leaving $a_{jj}(t)\, w(t)$ . Summing over $j$ ,

w'(t) = \bigl(\operatorname{tr}A(t)\bigr)\, w(t), \qquad w(0) = 1 ,

a scalar linear equation whose solution is $w(t) = \exp\int_0^t \operatorname{tr}A$ . At $t = T$ this is the stated formula, and since an exponential never vanishes, $\det M \neq 0$ .

Step 3: independence of the fundamental matrix. By Linear systems and the matrix exponential, any fundamental matrix is $\Psi(t) = \Phi(t)Q$ for the constant nonsingular $Q = \Psi(0)$ . Then $\Psi(t+T) = \Phi(t+T)Q = \Phi(t)MQ = \Psi(t)\,Q^{-1}MQ$ , so $M_\Psi = Q^{-1}MQ$ is similar to $M$ and has the same eigenvalues with the same multiplicities.

Step 4: a logarithm of $M$ . We claim that every nonsingular $M \in \mathbb{C}^{n\times n}$ is $e^{B}$ for some $B$ . Put $M$ in Jordan form (Jordan normal form), $M = SJS^{-1}$ with $J$ the direct sum of blocks $J_k = \lambda_k(I + N_k)$ , where $N_k$ is nilpotent and $\lambda_k \neq 0$ because $M$ is nonsingular. Fix any $\log\lambda_k$ with $e^{\log \lambda_k} = \lambda_k$ and set

L_k = (\log \lambda_k)\, I + \ell(N_k), \qquad \ell(X) = \sum_{m \ge 1} \frac{(-1)^{m+1}}{m}\, X^m ,

the sum in $\ell$ being finite because $N_k$ is nilpotent. In the ring of formal power series the identity $\exp(\ell(X)) = 1 + X$ holds, and both sides terminate when $X$ is replaced by a nilpotent matrix, so $e^{\ell(N_k)} = I + N_k$ . The two summands in $L_k$ commute (one is a scalar matrix), so $e^{L_k} = \lambda_k (I + N_k) = J_k$ by the product rule for commuting exponentials (Functions of matrices and the matrix exponential). Taking $B = S\bigl(\bigoplus_k L_k\bigr)S^{-1}$ and using $e^{SXS^{-1}} = Se^{X}S^{-1}$ gives $e^{B} = M$ .

Step 5: the normal form. Set $R = B/T$ , so $e^{TR} = M$ , and define $P(t) = \Phi(t)e^{-tR}$ . Then $P$ is continuously differentiable and nonsingular (a product of nonsingular matrices), $P(0) = I_n$ , and

P(t + T) = \Phi(t+T)\,e^{-(t+T)R} = \Phi(t)\,M\,e^{-TR}e^{-tR} = \Phi(t)e^{-tR} = P(t) ,

using Step 1, then $e^{-(t+T)R} = e^{-TR}e^{-tR}$ (the exponents commute), then $Me^{-TR} = I_n$ . So $P$ has period $T$ and $\Phi(t) = P(t)e^{tR}$ . For the change of variable, let $x$ solve the system and put $y = P(t)^{-1}x$ ; since $x(t) = \Phi(t)c = P(t)e^{tR}c$ for a constant $c$ , we get $y(t) = e^{tR}c$ , which satisfies $y' = Ry$ . Conversely $x = P(t)y$ recovers a solution. Finally, the eigenvalues of $R$ are exponents of the multipliers, with matching multiplicities. Indeed $\ell(N_k)$ is a sum of positive powers of a nilpotent matrix and so is nilpotent, hence $L_k = (\log\lambda_k)I + \ell(N_k)$ has $\log\lambda_k$ as its only eigenvalue, with multiplicity the size of the block $J_k$ . Therefore the eigenvalues of $R = S\bigl(\bigoplus_k L_k/T\bigr)S^{-1}$ , listed with algebraic multiplicity, are the numbers $(\log\lambda_k)/T$ repeated according to block size, matching the multipliers $\lambda_k$ repeated the same way; and $e^{T \cdot (\log\lambda_k)/T} = \lambda_k$ .

Step 6: Floquet solutions. Let $Mv = \rho v$ with $v \neq 0$ and put $x(t) = \Phi(t)v$ . This is a non-trivial solution, and $x(t+T) = \Phi(t+T)v = \Phi(t)Mv = \rho\, x(t)$ by Step 1. Choose $\mu$ with $e^{\mu T} = \rho$ and set $p(t) = e^{-\mu t}x(t)$ ; then $p(t+T) = e^{-\mu t}e^{-\mu T}\rho\, x(t) = p(t)$ , so $p$ has period $T$ and $x(t) = e^{\mu t}p(t)$ . Conversely any $x$ of that form satisfies $x(t+T) = \rho x(t)$ . This proves part 5 and completes Floquet's theorem. $\square$

Proof(What the multipliers decide)Original

Write any $t \ge 0$ as $t = kT + s$ with $k \in \mathbb{N} \cup \{0\}$ and $s \in [0, T)$ . A solution is $x(t) = \Phi(t)c$ with $c = x(0)$ , so by Step 1 above

x(kT + s) = \Phi(s)\, M^{k} c .

On the compact interval $[0,T]$ the continuous functions $s \mapsto \|\Phi(s)\|$ and $s \mapsto \|\Phi(s)^{-1}\|$ attain finite maxima $K$ and $K'$ , and $K' > 0$ . Hence

\frac{1}{K'}\,\bigl\|M^{k}c\bigr\| \;\le\; \|x(kT+s)\| \;\le\; K\,\bigl\|M^{k}c\bigr\| ,

the lower bound because $\|M^kc\| = \|\Phi(s)^{-1}x(kT+s)\| \le K'\|x(kT+s)\|$ . So the behaviour of all solutions as $t \to +\infty$ is the behaviour of the powers $M^k$ , and nothing else.

Powers of $M$ . Use the Jordan form again: $M^k = S J^k S^{-1}$ , and a block $\lambda(I+N)$ of size $d$ has

\bigl(\lambda(I+N)\bigr)^{k} = \lambda^{k}\sum_{m=0}^{d-1}\binom{k}{m} N^{m} ,

a finite sum because $N^d = 0$ . Writing it out, the entry of the block in the $m$ -th superdiagonal is $\lambda^{k-m}\binom{k}{m}$ for $0 \le m \le d-1$ , and everything below the diagonal is zero. For fixed $m$ the binomial coefficient $\binom{k}{m} = k(k-1)\cdots(k-m+1)/m!$ is a polynomial in $k$ of degree exactly $m$ , so each such entry has modulus comparable to $|\lambda|^{k}k^{m}$ for large $k$ . Three cases follow. If $|\lambda| < 1$ every entry of the block tends to $0$ , since a geometric factor beats a polynomial one. If $|\lambda| = 1$ the diagonal entries have modulus $1$ and do not tend to $0$ , while the entry at $m = 1$ , present as soon as $d \ge 2$ , has modulus $k$ and is unbounded. If $|\lambda| > 1$ the diagonal entries already blow up. Collecting the blocks: $J^{k} \to 0$ exactly when every eigenvalue has $|\lambda| < 1$ , and $\{J^{k}\}$ is bounded exactly when every eigenvalue has $|\lambda| < 1$ or has $|\lambda| = 1$ with all of its blocks of size $d = 1$ — that is, when the modulus-one eigenvalues are semisimple. Since $M^{k} = SJ^{k}S^{-1}$ with $S$ fixed and invertible, the same holds for $M^{k}$ , and since $\mathbb{C}^n$ has a basis, $\{M^{k}c\}$ bounded for every $c$ is the same as $\{M^{k}\}$ bounded. Parts 1, 2 and 3 follow, part 3 because an eigenvector $v$ for $|\rho| > 1$ gives $\|M^kv\| = |\rho|^k\|v\| \to \infty$ . For the rate in part 1, choose $r$ with $\max_j|\rho_j| < r < 1$ ; then $\|M^k\| \le C_0 r^{k}$ for some $C_0$ by the display above, and $\|x(kT+s)\| \le KC_0 r^{k}\|c\| \le (KC_0/r)\, r^{\,t/T}\|c\|$ since $k \ge t/T - 1$ .

Part 4. If $x$ is a non-trivial solution of period $T$ then $x(0) = x(T) = Mx(0)$ , and $x(0) \neq 0$ by uniqueness, so $1$ is an eigenvalue of $M$ . Conversely an eigenvector for $\rho = 1$ gives a Floquet solution with $x(t+T) = x(t)$ by Step 6. This is the matrix form of the criterion that a solution is periodic exactly when its value and derivative repeat after one period (Coddington (1989), Ch. 3, §3, Exercise 6(b), p. 110). The same argument with $x(T) = -x(0)$ handles $\rho = -1$ .

Part 5. For $y'' + \alpha(t)y = 0$ the coefficient matrix has zero trace, so $\det M = e^{0} = 1$ by Step 2; equivalently the Wronskian of the normalised basis is identically $1$ (Coddington (1989), Ch. 3, §4, Exercise 3(a), p. 116). A $2\times2$ matrix satisfies $\rho^2 - (\operatorname{tr}M)\rho + \det M = 0$ , which here is $\rho^2 - \tau\rho + 1 = 0$ , so $\rho_1\rho_2 = 1$ .

If $|\tau| < 2$ the discriminant $\tau^2 - 4$ is negative, so the roots are a conjugate pair of modulus $\sqrt{\rho_1\rho_2} = 1$ , and they are distinct, hence semisimple. Part 2 gives boundedness on $[0,\infty)$ ; applying the same argument to the time-reversed system, or using $\Phi(t - kT) = \Phi(t)M^{-k}$ with $M^{-1}$ having the same two multipliers inverted and still of modulus $1$ , gives boundedness on $(-\infty, 0]$ as well.
If $|\tau| > 2$ the roots are real, distinct and reciprocal, so one has modulus $> 1$ : part 3 applies.
$\rho = 1$ is a root exactly when $1 - \tau + 1 = 0$ , that is $\tau = 2$ , and $\rho = -1$ exactly when $1 + \tau + 1 = 0$ , that is $\tau = -2$ ; by part 4 these are the conditions for a non-trivial periodic, respectively antiperiodic, solution. In the normalised basis $\tau = \varphi_1(T) + \varphi_2'(T)$ , because that basis is exactly the one whose matrix is $\Phi$ . At $\tau = 2$ the only multiplier is $1$ with algebraic multiplicity $2$ , so by part 2 all solutions are bounded precisely when $M$ is semisimple, that is $M = I_2$ ; similarly for $\tau = -2$ and $M = -I_2$ . $\square$

Remark(Where the hypotheses are used, and the trap that catches everyone)Original

Continuity of $A$ on all of $\mathbb{R}$ is what makes $\Phi$ exist on the whole line, so that $\Phi(T)$ is defined and $\Phi(t+T)$ makes sense (Existence and uniqueness of solutions). Only continuity is needed, not smoothness: nothing in the proof differentiates $A$ .
Periodicity with the same $T$ in every entry enters exactly once, in Step 1. It is indispensable, not technical: if $A(t) = A_1(t) + A_2(t)$ with $A_1$ of period $1$ and $A_2$ of period $\sqrt 2$ , there is no $T$ at all, no monodromy matrix, and no multipliers. Quasi-periodic systems are a genuinely harder subject.
Linearity is used in Step 3 (every fundamental matrix is $\Phi Q$ ) and throughout the second proof, where the solution is $\Phi(t)c$ with $c$ constant.
Invertibility of $M$ comes from Step 2 and is needed for the multipliers to have logarithms in Step 4. Without it the normal form is false: $e^{TR}$ is always invertible.
$T$ is a choice. Replacing $T$ by $2T$ replaces every multiplier $\rho$ by $\rho^2$ , so a system with multipliers $\pm 1$ at period $T$ has both multipliers equal to $1$ at period $2T$ . The decay and boundedness criteria are unaffected, since $|\rho| < 1$ if and only if $|\rho^2| < 1$ , but the periodicity criterion in part 4 is not: it detects solutions of the period you tested with.
The trap: the eigenvalues of $A(t)$ tell you nothing. It is tempting to freeze $t$ , compute the eigenvalues of $A(t)$ and conclude stability when they all have negative real part for every $t$ . This is false. Take $A(t) = \begin{pmatrix} -1 + \tfrac32\cos^2 t & 1 - \tfrac32 \sin t\cos t \\[2pt] -1 - \tfrac32 \sin t \cos t & -1 + \tfrac32 \sin^2 t \end{pmatrix} ,$ which has period $T = \pi$ . Its trace is $-\tfrac12$ and its determinant is $\tfrac12$ at every $t$ , so its eigenvalues are the constants $(-1 \pm i\sqrt7)/4$ , with real part $-1/4 < 0$ for all $t$ . Yet $x_1(t) = e^{t/2}\begin{pmatrix} \cos t \\ -\sin t\end{pmatrix}, \qquad x_2(t) = e^{-t}\begin{pmatrix} \sin t \\ \cos t \end{pmatrix}$ are solutions, as substitution confirms, and they are the columns of the principal fundamental matrix because $x_1(0) = e_1$ and $x_2(0) = e_2$ . Hence $M = \Phi(\pi) = \begin{pmatrix} -e^{\pi/2} & 0 \\ 0 & -e^{-\pi}\end{pmatrix} ,$ with multipliers $-e^{\pi/2} \approx -4.8105$ and $-e^{-\pi} \approx -0.0432$ . The first has modulus greater than $1$ , so every solution whose initial vector is not a multiple of $e_2$ blows up. As a check, $\det M = e^{-\pi/2}$ agrees with $\exp\int_0^\pi \operatorname{tr}A = e^{-\pi/2}$ from Step 2. The instantaneous eigenvalues describe only the frozen system $x' = A(t_0)x$ . The true solution must in addition follow eigenvectors that rotate as $t$ advances — here those of $A(t)$ turn once per period — and that rotation feeds in growth the frozen picture cannot see. Only the monodromy matrix, which integrates over a whole period, accounts for both effects.

Remark(Real coefficients and the doubled period)Original

Suppose $A$ is real. Then $\Phi$ and $M$ are real, but $R$ in part 4 of Floquet's theorem need not be, and no amount of cleverness can make it so. Note first that any normal form $\Phi(t) = P(t)e^{tR}$ of period $T$ has $P(0) = \Phi(0) = I_n$ and therefore $M = \Phi(T) = P(T)e^{TR} = P(0)e^{TR} = e^{TR}$ , so a real normal form of period $T$ exists only if $M$ is the exponential of a real matrix. In the example just given, $M = \operatorname{diag}(-e^{\pi/2}, -e^{-\pi})$ . If $M$ were $e^{B}$ with $B = TR$ a real $2\times 2$ matrix, the eigenvalues $\beta_1, \beta_2$ of $B$ would satisfy $e^{\beta_i} < 0$ and so could not be real; a real $2 \times 2$ matrix with no real eigenvalue has a conjugate pair, and then $e^{\beta_2} = \overline{e^{\beta_1}}$ would force the two eigenvalues of $M$ to be complex conjugates — but they are distinct real numbers. So a real normal form of period $T$ does not exist here.

Doubling the period repairs this: there are a real $2T$ -periodic nonsingular $P$ and a real constant $R$ with $\Phi(t) = P(t)e^{tR}$ . It suffices to produce a real $S$ with $e^{S} = M^2$ , for then $R = S/(2T)$ and $P(t) = \Phi(t)e^{-tR}$ are real, and, exactly as in Step 5 but with $\Phi(t+2T) = \Phi(t)M^2$ and $e^{2TR} = M^2$ , we get $P(t+2T) = \Phi(t+2T)e^{-(t+2T)R} = \Phi(t)M^2e^{-2TR}e^{-tR} = P(t)$ . To build $S$ , split $\mathbb{C}^n$ into the generalised eigenspaces $E_\lambda$ of $M$ and write $M = \lambda(I + N_\lambda)$ on $E_\lambda$ with $N_\lambda$ nilpotent, using the series $\ell$ of Step 4.

If $\lambda > 0$ , set $S = 2(\log\lambda)I + 2\ell(N_\lambda)$ on $E_\lambda$ , with the real logarithm; then $e^{S} = \lambda^2(I+N_\lambda)^2 = M^2$ there.
If $\lambda < 0$ , note that $M^2 = \lambda^2(I + N'_\lambda)$ on $E_\lambda$ with $N'_\lambda = 2N_\lambda + N_\lambda^2$ nilpotent, and set $S = (\log \lambda^2)I + \ell(N'_\lambda)$ , again with a real logarithm. This is where the squaring is needed: a negative $\lambda$ has no real logarithm, but $\lambda^2$ does.
If $\lambda \notin \mathbb{R}$ , pick one eigenvalue out of each conjugate pair, call it $\lambda$ , fix any $\log\lambda$ and set $S = 2(\log\lambda)I + 2\ell(N_\lambda)$ on $E_\lambda$ ; on $E_{\bar\lambda}$ define $S$ by $Sv = \overline{S\bar v}$ , which is legitimate because complex conjugation maps $E_{\bar\lambda}$ onto $E_\lambda$ ( $M$ being real), and which gives $e^{S}v = \overline{M^2\bar v} = M^2 v$ there.

For real $\lambda$ the subspace $E_\lambda = \ker (M - \lambda I)^n$ is carried to itself by conjugation and the formulas above are real, so in every case $\overline{Sv} = S\bar v$ for all $v$ . A matrix commuting with conjugation in this way has real entries, and $e^{S} = M^2$ block by block. The price of realness is that $P$ repeats only after $2T$ ; a real Floquet solution attached to a negative multiplier $\rho$ is $e^{\sigma t}q(t)$ with $\sigma = T^{-1}\log|\rho|$ and $q(t+T) = -q(t)$ , so it has the doubled period built in.

Method(Finding the multipliers and reading off the stability)Original

Check the hypotheses. Confirm that $A$ is continuous and that a single $T > 0$ satisfies $A(t+T) = A(t)$ entrywise. For a scalar equation $y^{(n)} + a_1(t)y^{(n-1)} + \cdots + a_n(t)y = 0$ with $T$ -periodic coefficients, first rewrite it as a system (Linear systems and the matrix exponential); the companion matrix inherits the period.
Integrate one period, $n$ times. Solve $x' = A(t)x$ on $[0, T]$ with each standard basis vector $e_1, \dots, e_n$ as initial data, exactly if the system allows it and numerically otherwise. Assemble the results as the columns of $M$ . One period of integration is all that is ever required.
Check the determinant. Compare $\det M$ with $\exp\int_0^T \operatorname{tr}A$ . For a numerical integration this is a cheap and sharp error estimate; for Hill's equation it is the test $\det M = 1$ .
Compute the multipliers $\rho_j$ as the eigenvalues of $M$ , and the exponents as $\mu_j = T^{-1}\log \rho_j$ , remembering that each $\mu_j$ is fixed only up to adding a multiple of $2\pi i/T$ and that only $\operatorname{Re}\mu_j = T^{-1}\log|\rho_j|$ is unambiguous.
Read the verdict. All $|\rho_j| < 1$ : every solution decays exponentially. Some $|\rho_j| > 1$ : there are exponentially growing solutions. All $|\rho_j| \le 1$ with the modulus-one ones semisimple: all solutions bounded, none decaying to zero in general. A multiplier equal to $1$ or $-1$ : a periodic or antiperiodic solution exists, and you are on a stability boundary.
In dimension two with $\operatorname{tr}A \equiv 0$ , skip the eigenvalue computation: $\det M = 1$ , so the single number $\tau = \operatorname{tr}M$ decides everything by part 5 of the second statement. The stability diagram of a parametrised family is the set where $|\tau| < 2$ , and its boundary is the curve $|\tau| = 2$ .
To get a Floquet solution, take an eigenvector $v$ of $M$ for the multiplier of interest and solve the system once more with $x(0) = v$ ; the resulting solution reproduces itself multiplied by $\rho$ after each period.

Example(A square-wave Hill equation)Original

Take the equation

y'' + \alpha(t)\,y = 0, \qquad \alpha(t) = \begin{cases} 4, & 0 \le t < 1, \\ 1, & 1 \le t < 2,\end{cases} \qquad \alpha(t+2) = \alpha(t) ,

a stiffness that is switched between two values twice per period $T = 2$ . The coefficient is piecewise constant rather than continuous, so by a solution we mean a $C^1$ function satisfying the equation on each interval of constancy. Nothing in the theory above changes: the initial-value problem still has a unique solution, obtained by solving on $[0,1)$ and using the value and derivative at $t = 1$ as new initial data, the solution still depends linearly on $(y(0), y'(0))$ , and Step 1 only used $A(t+T) = A(t)$ together with uniqueness.

The monodromy matrix. On an interval of length $h$ where $\alpha \equiv \omega^2$ with $\omega > 0$ , the general solution is a combination of $\cos \omega t$ and $\sin\omega t$ , so the matrix carrying $(y, y')$ across that interval is

C(\omega, h) = \begin{pmatrix} \cos\omega h & \omega^{-1}\sin\omega h \\ -\omega \sin\omega h & \cos\omega h\end{pmatrix} ,

whose determinant is $\cos^2 + \sin^2 = 1$ . Our period consists of $h = 1$ at $\omega_1 = 2$ followed by $h = 1$ at $\omega_2 = 1$ , and the second step acts after the first, so

M = C(1, 1)\, C(2, 1) .

Multiplying out, with $c_j = \cos \omega_j$ and $s_j = \sin\omega_j$ ,

M = \begin{pmatrix} c_1 c_2 - 2 s_1 s_2 & \tfrac12 s_1 c_2 + s_2 c_1 \\[2pt] -s_2 c_1 - 2 s_1 c_2 & c_1c_2 - \tfrac12 s_1 s_2 \end{pmatrix} = \begin{pmatrix} -1.755140 & -0.104528 \\ -0.632416 & -0.607419 \end{pmatrix} ,

where $\cos 1 = 0.540302$ , $\sin 1 = 0.841471$ , $\cos 2 = -0.416147$ , $\sin 2 = 0.909297$ . Its determinant is $1$ , as it must be, since $\det C = 1$ for each factor and $\operatorname{tr}A \equiv 0$ .

The trace test. The trace is

\tau = 2\cos 1 \cos 2 - \Bigl(\frac{\omega_1}{\omega_2} + \frac{\omega_2}{\omega_1}\Bigr)\sin 1 \sin 2 = 2\cos 1\cos 2 - \tfrac52 \sin 1 \sin 2 = -2.362559 .

Since $|\tau| > 2$ , part 5 of the second statement says the equation has exponentially growing solutions: the switching destabilises an oscillator that is stable at either fixed stiffness. The multipliers are the roots of $\rho^2 + 2.362559\rho + 1 = 0$ ,

\rho_{1} = \frac{\tau - \sqrt{\tau^2-4}}{2} = -1.810104, \qquad \rho_2 = \frac{\tau+\sqrt{\tau^2-4}}{2} = -0.552454 ,

and indeed $\rho_1\rho_2 = 1.000000$ . The growth rate is

\sigma = \frac{\log|\rho_1|}{T} = \frac{\log 1.810104}{2} = 0.296692 ,

so the unstable Floquet solution has amplitude $e^{\sigma t}$ , multiplying by $1.810104$ every $2$ units of time. Because $\rho_1 < 0$ , that solution is $e^{\sigma t}q(t)$ with $q(t+2) = -q(t)$ , of period $4$ .

Checking the number. An eigenvector for $\rho_1$ comes from the first row of $M - \rho_1 I$ : taking $v_1 = 1$ ,

v_2 = \frac{\rho_1 - M_{11}}{M_{12}} = \frac{-1.810104 + 1.755140}{-0.104528} = 0.525836 .

Start the equation at $y(0) = 1$ , $y'(0) = 0.525836$ . On $[0,1]$ , $y = \cos 2t + 0.262918\sin 2t$ , so $y(1) = \cos 2 + 0.262918\sin 2 = -0.416147 + 0.239071 = -0.177076$ and $y'(1) = -2\sin 2 + 0.525836\cos 2 = -1.818595 - 0.218825 = -2.037420$ . On $[1,2]$ the frequency is $1$ , so $y(t) = y(1)\cos(t-1) + y'(1)\sin(t-1)$ and

y(2) = -0.177076\,(0.540302) + (-2.037420)(0.841471) = -0.095675 - 1.714430 = -1.810105 .

The predicted value is $\rho_1 y(0) = -1.810104$ , and the last digit of the computed $-1.810105$ is the rounding carried by the six-digit intermediate values. The derivative comes out the same way: $y'(2) = -y(1)\sin 1 + y'(1)\cos 1 = -0.951818$ , against the predicted $\rho_1 v_2 = -1.810104\,(0.525836) = -0.951818$ .

A stable neighbour. Halving both intervals, $\alpha = 4$ on $[0,\tfrac12)$ and $\alpha = 1$ on $[\tfrac12, 1)$ with $T = 1$ , the same formula gives $\tau = 2\cos\tfrac12\cos 1 - \tfrac52 \sin\tfrac12 \sin 1 = -0.060237$ . Now $|\tau| < 2$ , the multipliers are $-0.030118 \pm 0.999546\,i$ of modulus exactly $1$ , and every solution is bounded. Switching the same two stiffnesses twice as fast turns instability into stability — which is the whole content of a parametric stability chart.

ApplicationOriginal

Parametric resonance. A pendulum whose support is driven vertically, or a child pumping a swing by changing the effective length twice per swing, obeys Hill's equation $y'' + \alpha(t)y = 0$ with $\alpha$ periodic; the smooth model is Mathieu's equation $y'' + (\delta + \varepsilon\cos 2t)y = 0$ . The theory above says that the $(\delta, \varepsilon)$ plane splits into regions $|\tau(\delta,\varepsilon)| < 2$ where the swing stays small and regions $|\tau| > 2$ where it grows exponentially, separated by the curves $\tau = \pm 2$ on which a $T$ -periodic or $2T$ -periodic solution exists. Those curves are the familiar tongues issuing from $\delta = m^2$ , and the tongue at $\delta = 1$ — pumping at twice the natural frequency — is the one the child uses. The same equation with the same analysis governs the radio-frequency quadrupole (Paul) ion trap, where the stable region of the chart is precisely the set of ion masses the trap confines.
Stability of limit cycles. If $\gamma$ is a $T$ -periodic solution of an autonomous nonlinear system $x' = f(x)$ , the equation for a small deviation is the $T$ -periodic linear system $v' = Df(\gamma(t))\,v$ , so its Floquet multipliers measure the growth of perturbations. One multiplier is always $1$ : differentiating $\gamma'(t) = f(\gamma(t))$ shows that $v = \gamma'$ solves the variational equation, and $\gamma'$ has period $T$ . Strogatz (Strogatz (1994), §8.7, p. 282) calls this the trivial multiplier, discards it as mere time translation, and states that the closed orbit is linearly stable exactly when the remaining $n-1$ multipliers of the linearized Poincaré map have modulus less than $1$ — which is part 1 of the statement above, read on the complement of the direction along the orbit.
Electrons in a one-dimensional crystal. The stationary Schrödinger equation $-\psi'' + V(t)\psi = E\psi$ with $V$ of period $T$ is Hill's equation with $\alpha = E - V$ . For each energy $E$ the trace $\tau(E)$ of the monodromy matrix is a continuous function, and no solution grows exponentially exactly when $|\tau(E)| \le 2$ . The energy intervals where this holds are the allowed bands of the solid, the intervals where $|\tau(E)| > 2$ are the forbidden gaps, and the band edges are the energies with $\tau(E) = \pm 2$ , where the Floquet solutions are genuinely periodic or antiperiodic. The Floquet solution $e^{\mu t}p(t)$ with $p$ periodic is the Bloch wave, and $\mu = ik$ defines the crystal momentum inside a band.
Beam optics. In an alternating-gradient accelerator the transverse displacement of a particle satisfies Hill's equation with a focusing strength that repeats each lattice cell. Designers compute the one-cell transfer matrix — our $M$ — and require $|\operatorname{tr}M| < 2$ , writing $\operatorname{tr}M = 2\cos\Theta$ to read off the phase advance $\Theta$ per cell. The stability criterion is exactly part 5 above.

Method

Solving initial-value problems with the Laplace transform

Adapted

Add to PDF

Turns a constant-coefficient initial-value problem into one algebraic equation for the transform of the unknown, the subsidiary equation p(s)Y(s) = Q(s) + F(s), whose numerator carries the initial data, and shows how partial fractions, the shifting theorem and the convolution rule recover the solution, including for forcing that is switched on at a later time or delivered as an impulse.

Prerequisites: The Laplace transform, Constant-coefficient equations, Partial fractions for inverse transforms, The Heaviside step and related functions, Delta-function identities, Test functions and distributions, Green's function for a boundary-value problem

Note(Notation, and the three rules the method uses)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.8, Eqs. (20.147)-(20.149), p. 1016 and Eq. (20.166), p. 1023; Riley, Hobson & Bence (2006), §13.2, Eq. (13.53), p. 453 and §13.2.1, Eqs. (13.57)-(13.58), pp. 455-456; Strauss (2008), §12.5, Eq. (1), p. 353 and the property table (ii), (iii), (vii), (ix), p. 354

Throughout this unit a function of $t$ carries a lower-case letter and its transform the corresponding capital:

\mathcal{L}\{v\}(s) = V(s) = \int_0^\infty e^{-st}\, v(t)\,dt

(the Laplace transform); so the unknown $y$ has transform $Y$ , the forcing $f$ has transform $F$ , and so on.

The cited books do not agree on this, and the disagreement is pure notation. Strauss writes exactly the convention above, $F(s)$ for the transform of $f(t)$ ((Strauss (2008), §12.5, Eq. (1), p. 353)). Riley, Hobson and Bence keep the letter and put a bar on it, $\bar f(s)$ ((Riley, Hobson & Bence (2006), §13.2, Eq. (13.53), p. 453)). Arfken uses the opposite case to ours: for him the time function is the capital and the transform the lower-case letter, $\mathcal{L}\{F(t)\} = f(s)$ ((Arfken, Weber & Harris (2013), §20.8, Eqs. (20.147)-(20.149), p. 1016)). Every formula imported from Arfken below is therefore restated in the convention above, his $F$ becoming our $v$ and his $f$ our $V$ .

Three properties of the transform, proved in the unit cited above, are the whole machinery of this unit.

Derivatives. Integration by parts gives

\mathcal{L}\{v'\}(s) = s\,V(s) - v(+0) ,

and by iteration

\mathcal{L}\{v^{(n)}\}(s) = s^n V(s) - s^{n-1}v(+0) - s^{n-2}v'(+0) - \dots - v^{(n-1)}(+0) .

Arfken attaches three conditions to these identities ((Arfken, Weber & Harris (2013), §20.8, the text following Eq. (20.147), p. 1016)). The value at the origin is to be read as the one-sided limit $v(+0)$ . The derivative is allowed jumps, being asked only to be piecewise continuous on the whole half-line $t \ge 0$ . And the transform integrals of $v$ and of $v'$ must both converge rather than run off to infinity. Riley, Hobson and Bence derive the same identities and append the proviso $s > 0$ ((Riley, Hobson & Bence (2006), §13.2.1, Eqs. (13.57)-(13.58), pp. 455-456)); the remark on limitations below returns to the difference. 2. Shifting in $t$ (the Heaviside shifting theorem). For $b > 0$ ,

\mathcal{L}\bigl\{ H(t-b)\,v(t-b)\bigr\}(s) = e^{-bs}\,V(s) ,

with $H$ the Heaviside step. Arfken derives it in the form $\mathcal{L}\{v(t-b)\} = e^{-bs}V(s)$ on the convention that $v$ vanishes for negative argument, and notes that it is a Heaviside factor inserted into the integrand that enforces that convention; he calls the result the Heaviside shifting theorem ((Arfken, Weber & Harris (2013), §20.8, Eqs. (20.164)-(20.166), pp. 1022-1023)). Strauss lists it in the form above ((Strauss (2008), §12.5, property (vii), p. 354)). 3. Convolution. $\mathcal{L}\bigl\{\int_0^t v(t-\tau)\,w(\tau)\,d\tau\bigr\}(s) = V(s)\,W(s)$ ((Strauss (2008), §12.5, property (ix), p. 354)).

Two transform pairs are used constantly: $\mathcal{L}\{H(t-b)\} = e^{-bs}/s$ and $\mathcal{L}\{\delta(t-b)\} = e^{-bs}$ , both tabulated by Strauss and by Riley, Hobson and Bence ((Strauss (2008), §12.5, entries (8) and (9), p. 353); (Riley, Hobson & Bence (2006), §13.2, Table 13.1, p. 455)). The second is a statement about the delta distribution rather than about a function, and is read as $\int_0^\infty e^{-st}\delta(t-b)\,dt = e^{-bs}$ .

Method(The subsidiary equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.1.5, pp. 501-503; Arfken, Weber & Harris (2013), §20.8, pp. 1016-1017

Let $a_1, \dots, a_n$ be constants, let

p(s) = s^n + a_1 s^{n-1} + \dots + a_{n-1}s + a_n

be the characteristic polynomial of the operator $L = \frac{d^n}{dt^n} + a_1\frac{d^{n-1}}{dt^{n-1}} + \dots + a_n$ , and consider

L y = f(t), \qquad y(0) = y_0, \ y'(0) = y_1, \ \dots, \ y^{(n-1)}(0) = y_{n-1} .

Assume $y$ and its first $n$ derivatives have Laplace transforms in the sense of rule 1 above, and write, in the capital-for-transform convention of the note, $Y = \mathcal{L}\{y\}$ and $F = \mathcal{L}\{f\}$ . Then $Y$ satisfies one algebraic equation, the subsidiary equation

p(s)\,Y(s) = Q(s) + F(s), \qquad Q(s) = \sum_{m=1}^{n} a_{n-m}\sum_{j=0}^{m-1} s^{\,m-1-j}\, y_j \quad (a_0 := 1) ,

so that

Y(s) = \frac{Q(s)}{p(s)} + \frac{F(s)}{p(s)} .

Here $Q$ is a polynomial of degree at most $n-1$ built from the initial data alone, and $1/p(s)$ — the transfer function — knows only the equation. Consequently:

Separation of the two jobs. $y = y_{\text{free}} + y_{\text{forced}}$ with $Y_{\text{free}} = Q/p$ and $Y_{\text{forced}} = F/p$ . The first depends only on the initial data, the second only on the forcing.
The free part is the complementary function. $Q/p$ is a proper rational function, so partial fractions (partial fractions for inverse transforms) express it as a sum of terms $c/(s-\lambda)^k$ with $\lambda$ a root of $p$ and $k$ at most its multiplicity; inverting term by term with $\mathcal{L}^{-1}\{(s-\lambda)^{-k}\} = t^{k-1}e^{\lambda t}/(k-1)!$ gives exactly the solutions $t^{k-1}e^{\lambda t}$ of the constant-coefficient equation. The characteristic equation and the denominator of the transform are the same polynomial.
The forced part is a convolution. By rule 3, with $k = \mathcal{L}^{-1}\{1/p\}$ ,

y_{\text{forced}}(t) = \int_0^t k(t-\tau)\,f(\tau)\,d\tau .

The kernel $k$ is the impulse response: it is the solution of $Lk = 0$ with $k(0) = \dots = k^{(n-2)}(0) = 0$ , $k^{(n-1)}(0) = 1$ , and this integral is the Duhamel integral of nonhomogeneous systems in scalar form. Strauss's Example 2 is this statement for $n = 2$ : $u'' + \omega^2 u = f$ with zero data has $U(s) = F(s)/(s^2+\omega^2)$ , hence $u(t) = \int_0^t \omega^{-1}\sin\bigl(\omega(t-\tau)\bigr)f(\tau)\,d\tau$ ((Strauss (2008), §12.5, Example 2, pp. 354-355)).

For $n = 2$ the polynomial $Q$ is $Q(s) = (s + a_1)y_0 + y_1$ , which is the form Riley, Hobson and Bence produce in their worked example and the one met in practice ((Riley, Hobson & Bence (2006), §15.1.5, Eqs. (15.33)-(15.34), p. 502)).

The same manipulation handles a system of constant-coefficient equations: transforming each equation gives a set of simultaneous linear algebraic equations in the transforms, which are then solved by elimination, as in the coupled-circuit example of (Riley, Hobson & Bence (2006), §15.1.5, pp. 502-503) and the Earth-nutation example of (Arfken, Weber & Harris (2013), §20.8, Example 20.8.3, pp. 1018-1019). Arfken makes this the headline property of the method: coupled constant-coefficient differential equations transform into simultaneous linear algebraic equations ((Arfken, Weber & Harris (2013), §20.8, p. 1016)).

Idea of the proofOriginal

Each differentiation costs one factor of $s$ and leaves behind one piece of initial data; that is the whole of rule 1. Apply it to every term of a constant-coefficient equation and the differential operator becomes multiplication by its characteristic polynomial, while the debris of initial data collects into a polynomial of lower degree. One division then solves the problem in the transform variable, and all the remaining work — partial fractions, table look-up, shifting — is about getting back. The reason the initial data land in the numerator and the equation in the denominator is worth noticing: it is why changing the initial conditions cannot change which exponentials appear in the answer, only their coefficients.

ProofOriginal

Transform each term. Write $\mathcal{L}\{y^{(m)}\} = s^m Y - \sum_{j=0}^{m-1}s^{\,m-1-j}y_j$ , which is rule 1 with $n$ replaced by $m$ (and is an empty sum for $m=0$ , i.e. $\mathcal{L}\{y\} = Y$ ). The equation $Ly = f$ reads $\sum_{m=0}^{n} a_{n-m}\,y^{(m)} = f$ with $a_0 = 1$ , and $\mathcal{L}$ is linear, so

\sum_{m=0}^{n} a_{n-m}\Bigl[s^m Y - \sum_{j=0}^{m-1}s^{\,m-1-j}y_j\Bigr] = F(s) .

The coefficients being constants is what allows this: $\mathcal{L}$ is linear over constants, not over functions of $t$ , and that restriction is exactly the limitation discussed in the remark below.

Collect. $\sum_{m=0}^n a_{n-m}s^m = p(s)$ by definition, and the remaining double sum is $Q(s)$ , which therefore satisfies $p(s)Y(s) - Q(s) = F(s)$ . Each inner sum has highest power $s^{m-1}$ , and the outer index runs to $m = n$ , so $\deg Q \le n-1$ .
Part 1 is the linearity of $\mathcal{L}^{-1}$ applied to the splitting $Y = Q/p + F/p$ : the two inverse transforms are separately well defined, and their sum is $y$ because the inverse transform is linear and (up to the proviso recorded below) unique.
Part 2. $\deg Q \le n - 1 < n = \deg p$ , so $Q/p$ is proper and the partial-fraction theorem applies with no polynomial part. Each resulting term is $c\,(s-\lambda)^{-k}$ with $\lambda$ a root of $p$ and $1 \le k \le$ its multiplicity. From $\mathcal{L}\{t^{k-1}\} = (k-1)!/s^{k}$ and the substitution rule $\mathcal{L}\{e^{\lambda t}v(t)\}(s) = V(s - \lambda)$ ((Arfken, Weber & Harris (2013), §20.8, Eq. (20.157), p. 1020)) we get $\mathcal{L}\{t^{k-1}e^{\lambda t}\} = (k-1)!\,(s-\lambda)^{-k}$ , which inverts each term. The functions so produced are precisely a basis of the solution space of $Ly = 0$ .
Part 3. $1/p$ is proper, so it is the transform of some $k$ by step 4; then $F/p = \mathcal{L}\{k\}\mathcal{L}\{f\}$ , and rule 3 identifies the inverse transform as $\int_0^t k(t-\tau)f(\tau)d\tau$ . To see that $k$ has the stated initial data, apply the subsidiary equation to the problem $Lk = 0$ with those data: $Q(s)$ then reduces to its $m = n$ , $j = n-1$ term, $a_0 \cdot s^0 \cdot k^{(n-1)}(0) = 1$ , so $K(s) = 1/p(s)$ . Conversely a function with transform $1/p$ is unique, so $k$ is that solution. Finally $y_{\text{forced}}$ has $n-1$ vanishing derivatives at $0$ , because differentiating $\int_0^t k(t-\tau)f(\tau)d\tau$ $j < n$ times gives $\int_0^t k^{(j)}(t-\tau)f(\tau)\,d\tau$ plus terms containing $k^{(i)}(0)$ with $i < j \le n-1$ , all zero; so the splitting in part 1 matches the initial data as claimed. $\square$

Method(Solving a constant-coefficient initial-value problem by transform)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.1.5, the Solution method, p. 503

Transform the whole equation, using the derivative rule, and substitute the given values of $y(0), y'(0), \dots$ as you go. The initial conditions enter here and are never fitted afterwards — the reason Riley, Hobson and Bence describe the method as giving the solution in one step instead of requiring a complementary function and a particular integral to be found separately ((Riley, Hobson & Bence (2006), §15.1.5, p. 501)).
Solve the resulting algebraic equation for $Y(s)$ . For a system, solve the simultaneous algebraic equations by elimination. The result is always a rational function of $s$ , possibly multiplied by factors $e^{-bs}$ if the forcing is switched on at times $b > 0$ .
Put each rational factor into partial fractions (partial fractions). Factor the denominator first; a real quadratic factor with complex roots is best kept intact and handled by completing the square, $s^2 + 2\alpha s + \beta^2 = (s+\alpha)^2 + \omega^2$ with $\omega^2 = \beta^2 - \alpha^2$ , after which the numerator is arranged as $A(s+\alpha) + B\omega$ so that the two standard pairs $e^{-\alpha t}\cos\omega t$ and $e^{-\alpha t}\sin\omega t$ apply. Arfken's damped-oscillator example is this step done once ((Arfken, Weber & Harris (2013), §20.8, Example 20.8.5, p. 1021)).
Invert term by term from a table ((Riley, Hobson & Bence (2006), §13.2, Table 13.1, p. 455)). A factor $e^{-bs}$ is removed last, by the shifting theorem: if $\mathcal{L}^{-1}\{\Psi(s)\} = \psi(t)$ , then $\mathcal{L}^{-1}\{e^{-bs}\Psi(s)\} = H(t-b)\,\psi(t-b)$ . Where a factor of the transform has no convenient inverse, use the convolution rule instead and leave the answer as an integral.
If the table fails, two further devices often rescue the inversion: differentiating a known transform with respect to a parameter, since $\mathcal{L}\{(-t)^n v(t)\}(s) = V^{(n)}(s)$ — Arfken's rule, restated in our convention ((Arfken, Weber & Harris (2013), §20.8, Eqs. (20.173)-(20.174), pp. 1024-1025)) — and in the last resort the Bromwich inversion integral $v(t) = (2\pi i)^{-1}\int_{\alpha - i\infty}^{\alpha+i\infty}e^{st}V(s)\,ds$ ((Strauss (2008), §12.5, Eq. (13), p. 354)), evaluated by residues as in the Laplace transform.
Check. Verify the initial conditions directly on the answer; verify the differential equation at one convenient value of $t$ or symbolically; and check the limit $s \to \infty$ of the subsidiary equation against $y(0)$ , since $sY(s) \to y(+0)$ .

When to use it, and when not. The transform is the method of choice when the forcing is discontinuous, impulsive, or switched on at a later time, because a step or a delta transforms to something no more complicated than $e^{-bs}/s$ or $e^{-bs}$ , whereas handling them by undetermined coefficients requires splitting the time axis and matching at the junctions by hand. It is the wrong tool when the coefficients are not constant, when the data are given at two different points instead of all at $t = 0$ (a boundary-value problem — use Green's functions), or when the forcing has no elementary transform.

Method(Switched and impulsive forcing)Original

Let $L$ be as above with $n \ge 1$ , let $k = \mathcal{L}^{-1}\{1/p\}$ be the impulse response, and let $b > 0$ . For zero initial data:

Impulse. The problem $Ly = c\,\delta(t-b)$ has transform $Y(s) = c\,e^{-bs}/p(s)$ , hence

y(t) = c\,H(t-b)\,k(t-b) .

So $y$ vanishes until the impulse arrives and then evolves as a free solution; $y, y', \dots, y^{(n-2)}$ are continuous at $t = b$ , while $y^{(n-1)}$ jumps there by $c$ . Nothing but the jump distinguishes an impulse from a change of initial condition: an impulse of strength $c$ at time $b$ is the same as restarting the equation at $t = b$ with $y^{(n-1)}$ increased by $c$ . 2. Step. The problem $Ly = c\,H(t-b)$ has $Y(s) = c\,e^{-bs}/\bigl(s\,p(s)\bigr)$ , hence

y(t) = c\,H(t-b)\,\kappa(t-b), \qquad \kappa(t) = \int_0^t k(\tau)\,d\tau ,

so the step response $\kappa$ is the integral of the impulse response. Now $y, y', \dots, y^{(n-1)}$ are all continuous at $t = b$ and $y^{(n)}$ jumps by $c$ . 3. Several pieces. Because the problem is linear, a forcing assembled from steps, impulses and ordinary pieces is handled term by term and the results added; the shifting theorem supplies one $H(t - b_i)$ per switch-on time.

Part 1 is to be read distributionally: $\delta$ is not a function and $y$ is not a classical solution at $t = b$ (distributions). The content that survives in classical terms is exactly the jump statement, and that is how the formula is verified below. The remark after the proof records the $n = 2$ case as Arfken gives it.

Idea of the proofOriginal

The transform turns "wait until time $b$ , then act" into a multiplication by $e^{-bs}$ , and multiplication is something an algebraic expression can absorb without any case distinction. So all three parts are the single observation that the forcing's transform factors as $e^{-bs}$ times something already inverted. The classical content — which derivative jumps and by how much — is then read off by integrating the differential equation across $t = b$ : only the top derivative can absorb a delta, so only it jumps, and the amount is the strength of the delta.

ProofOriginal

Impulse. $\mathcal{L}\{c\,\delta(t-b)\} = c\,e^{-bs}$ , so the subsidiary equation with zero data gives $Y = c\,e^{-bs}/p(s) = c\,e^{-bs}\mathcal{L}\{k\}$ , and the shifting theorem (rule 2, available because $b>0$ ) inverts this to $c\,H(t-b)k(t-b)$ .

For the jump statement, argue classically. On $t < b$ the function is $0$ and on $t > b$ it is $c\,k(t-b)$ , which solves $Ly = 0$ there. Since $k(0) = \dots = k^{(n-2)}(0) = 0$ and $k^{(n-1)}(0) = 1$ , the one-sided limits at $b$ are

y^{(j)}(b^-) = 0 = y^{(j)}(b^+) \ \ (j \le n-2), \qquad y^{(n-1)}(b^-) = 0, \quad y^{(n-1)}(b^+) = c ,

which is continuity of the first $n-1$ functions and a jump of $c$ in the last. Conversely, integrating $Ly = c\,\delta(t-b)$ over $(b-\varepsilon, b+\varepsilon)$ : the term $y^{(n)}$ integrates to $y^{(n-1)}(b+\varepsilon) - y^{(n-1)}(b-\varepsilon)$ , every lower term integrates to something tending to $0$ with $\varepsilon$ because $y^{(j)}$ is bounded for $j \le n-1$ , and the right-hand side integrates to $c$ . So the jump in $y^{(n-1)}$ must be $c$ , and no lower derivative may jump, since a jump in $y^{(j)}$ with $j \le n-2$ would make $y^{(j+1)}$ contain a delta and $y^{(n)}$ a derivative of a delta, which the right-hand side does not have.

Step. $\mathcal{L}\{c\,H(t-b)\} = c\,e^{-bs}/s$ , so $Y = c\,e^{-bs}/(s\,p(s))$ . By rule 1 read backwards — or by the convolution rule with $\mathcal{L}\{1\} = 1/s$ — division of a transform by $s$ corresponds to integration from $0$ : $\mathcal{L}\{\int_0^t k\} = \mathcal{L}\{k\}/s$ , since $\kappa' = k$ and $\kappa(0) = 0$ give $\mathcal{L}\{k\} = s\mathcal{L}\{\kappa\} - 0$ . Hence $1/(s\,p(s)) = \mathcal{L}\{\kappa\}$ and the shifting theorem finishes it. The continuity statement follows from part 1 by integration: $\kappa^{(j)} = k^{(j-1)}$ for $j \ge 1$ , so $\kappa, \dots, \kappa^{(n-1)}$ all vanish at $0$ and $\kappa^{(n)}(0) = k^{(n-1)}(0) = 1$ .
Several pieces. $\mathcal{L}$ and $\mathcal{L}^{-1}$ are linear, and so is $L$ ; a finite sum of forcings therefore produces the sum of the responses. $\square$

Remark(The second-order case, as Arfken reads it)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.8, Example 20.8.4, pp. 1019-1020

Part 1 appears in the cited sources only in its second-order case, and Arfken reads it mechanically rather than algebraically. His operator is the bare $L = m\,d^2/dt^2$ and his forcing a multiple of the delta, so the single jump the statement predicts is a jump in the velocity: the impulse transfers its strength in units of linear momentum to the particle, instantaneously, and after that the particle coasts. He then runs the identical computation for a ballistic galvanometer, where the deposit is of angular momentum and its size is proportional to the charge $q$ passed ((Arfken, Weber & Harris (2013), §20.8, Example 20.8.4, pp. 1019-1020)). Translating to the statement above: $n = 2$ , the derivative that jumps is $y^{(n-1)} = y'$ , and the strength $c$ of the impulse is the transferred momentum divided by $m$ — with the one difference that Arfken delivers his kick at $t = 0$ , whereas the statement above keeps $b > 0$ so that the shifting theorem can be used. The general- $n$ form of the statement — one impulse, one jump, in the top derivative only — is the handbook's extension of that reading, not Arfken's.

Remark(Limitations, and the small print)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.8, pp. 1016, 1025-1026; Riley, Hobson & Bence (2006), §13.2, p. 454 and §13.2.1, pp. 455-456

Constant coefficients are essential. If a coefficient depends on $t$ the transform of the equation is no longer algebraic. Arfken's Example 20.8.7 shows what happens instead: Bessel's equation of order zero, written here as $t\,w'' + w' + t\,w = 0$ , transforms — using also $\mathcal{L}\{t\,v(t)\}(s) = -V'(s)$ — into the first-order differential equation $(s^2+1)W'(s) + s\,W(s) = 0$ for the transform $W = \mathcal{L}\{w\}$ , which happens to be solvable and yields $W(s) = C(s^2+1)^{-1/2}$ ((Arfken, Weber & Harris (2013), §20.8, Example 20.8.7, pp. 1025-1026); Arfken writes the same calculation with his $F$ for $w$ and his $f$ for $W$ ). He presents it as an interesting application of the transform-differentiation rule, and gives the general warning separately, where the derivative rules are set up, by pointing the reader to that very example for what may happen when the coefficients are not constant, and by locating the power of the transform in the collapse of a constant-coefficient equation to an algebraic one ((Arfken, Weber & Harris (2013), §20.8, p. 1016)). The further comment is ours and not his: the outcome here is a piece of luck — the transformed problem turned out easier than the original — and is not what to expect, since a variable coefficient in general trades one differential equation for another of no lower order.
The one-sided value $v(+0)$ , not $v(0)$ . The derivative rule carries the limit from the right, because the transform integral starts at $t = 0$ . For a problem whose forcing switches on at $t=0$ this matters: $y(+0)$ is the initial condition the transform sees, and a jump of the solution at the origin is invisible to it.
The region of convergence. Riley, Hobson and Bence state the derivative rules with the proviso $s > 0$ ((Riley, Hobson & Bence (2006), §13.2.1, Eqs. (13.57)-(13.58), pp. 455-456)), an inequality that holds for the bounded and polynomially growing functions transformed in that subsection but is not sufficient in general: for $v(t) = e^{2t}$ neither side of the first identity exists at $s = 1$ . Their own general discussion supplies the condition that is needed, namely that for each $v$ there is a number $s_0$ with the defining integral convergent for $s > s_0$ and divergent for $s \le s_0$ , and their table lists $s_0$ for every entry ((Riley, Hobson & Bence (2006), §13.2, p. 453 and Table 13.1, p. 455)). This handbook states the rule with Arfken's hypothesis — that the integrals do not diverge — and reads the inequality on $s$ off the table entry in each case. In practice no harm comes of it, because the final inversion is done by recognizing transforms and not by evaluating any integral at a particular $s$ .
Uniqueness of the inverse. Riley, Hobson and Bence say that the inverse transform may for every practical purpose be treated as unique, and attach the footnote that two functions which disagree at only finitely many isolated points share one transform ((Riley, Hobson & Bence (2006), §13.2, p. 454 and its footnote)). This is why the method is entitled to assert that a transform found in a table is the solution, and why the value of a switched solution exactly at the switching time is not determined by the transform.
What the method assumes about the answer. Step 1 of the proof presumes that $y$ and its derivatives have transforms. For a constant-coefficient problem this is harmless, since the solution exists by the constant-coefficient theory and is a combination of $t^je^{\lambda t}$ and of integrals of the forcing, all of exponential type if $f$ is; but the presumption is there, and strictly the transform method verifies a solution it has guessed rather than proving existence.

Example(A damped oscillator kicked twice: a switched step and an impulse)Original

Solve

y'' + 2y' + 5y = H(t-1) + 3\,\delta(t-2), \qquad y(0) = y'(0) = 0 .

Physically: an oscillator with natural frequency $\sqrt5$ and damping $2$ , at rest, has a constant force switched on at $t = 1$ and receives a sharp blow of strength $3$ at $t = 2$ .

1. Transform (step 1). With zero initial data the polynomial $Q$ vanishes, so $p(s) = s^2 + 2s + 5$ and

(s^2+2s+5)\,Y(s) = \frac{e^{-s}}{s} + 3e^{-2s} , \qquad Y(s) = e^{-s}\,\underbrace{\frac{1}{s(s^2+2s+5)}}_{= G(s)} + 3e^{-2s}\underbrace{\frac{1}{s^2+2s+5}}_{= K(s)} .

2. Partial fractions and completing the square (steps 2-3). Write $G(s) = A/s + (Bs+C)/(s^2+2s+5)$ . Multiplying by $s(s^2+2s+5)$ ,

1 = A(s^2+2s+5) + s(Bs + C) .

Setting $s = 0$ gives $A = 1/5$ . Comparing coefficients of $s^2$ gives $A + B = 0$ , so $B = -1/5$ ; comparing coefficients of $s$ gives $2A + C = 0$ , so $C = -2/5$ . Hence

G(s) = \frac15\left[\frac1s - \frac{s+2}{s^2+2s+5}\right] .

Completing the square, $s^2 + 2s + 5 = (s+1)^2 + 4$ , so $\omega = 2$ and the numerator is arranged as $s + 2 = (s+1) + \tfrac12 \cdot 2$ :

\frac{s+2}{(s+1)^2+4} = \frac{s+1}{(s+1)^2+4} + \frac12\cdot\frac{2}{(s+1)^2+4} .

3. Invert (step 4). From the standard pairs $e^{-\alpha t}\cos\omega t \leftrightarrow (s+\alpha)/((s+\alpha)^2+\omega^2)$ and $e^{-\alpha t}\sin\omega t \leftrightarrow \omega/((s+\alpha)^2+\omega^2)$ with $\alpha = 1$ , $\omega = 2$ ,

k(t) = \mathcal{L}^{-1}\{K\} = \tfrac12 e^{-t}\sin 2t , \qquad g(t) = \mathcal{L}^{-1}\{G\} = \tfrac15\Bigl[1 - e^{-t}\bigl(\cos 2t + \tfrac12 \sin 2t\bigr)\Bigr] .

Here $k$ is the impulse response and $g$ the step response, and the general statement $\kappa' = k$ is visible: differentiating $g$ ,

g'(t) = \tfrac15 e^{-t}\Bigl[\bigl(\cos 2t + \tfrac12\sin 2t\bigr) - \bigl(-2\sin 2t + \cos 2t\bigr)\Bigr] = \tfrac15 e^{-t}\cdot\tfrac52\sin 2t = k(t) .

Removing the exponential factors by the shifting theorem,

\boxed{\;y(t) = H(t-1)\,g(t-1) + 3\,H(t-2)\,k(t-2)\;}

that is, $y \equiv 0$ on $[0,1]$ ; $y(t) = g(t-1)$ on $[1,2]$ ; and $y(t) = g(t-1) + 3k(t-2)$ for $t \ge 2$ .

4. Checks (step 6). $y(0) = y'(0) = 0$ because nothing happens before $t=1$ . At the switch-on time $t = 1$ : $g(0) = \tfrac15[1 - 1] = 0$ and $g'(0) = k(0) = 0$ , so $y$ and $y'$ are continuous there and only $y''$ jumps, by $1$ — which is part 2 of the statement with $n = 2$ . At the impulse time $t = 2$ : $k(0) = 0$ but $k'(0) = \tfrac12(2) = 1$ , so $y$ is continuous while $y'$ jumps by $3 \cdot 1 = 3$ — part 1 of the statement. And as $t \to \infty$ all the exponentials die and $y \to \tfrac15$ , which is the static response $1/p(0) = 1/5$ to a constant force of $1$ , as it must be.

5. Numbers. At $t = 3$ , two units after the step and one after the blow,

y(3) = g(2) + 3k(1) = \tfrac15\bigl[1 - e^{-2}(\cos 4 + \tfrac12 \sin 4)\bigr] + \tfrac32 e^{-1}\sin 2 .

With $\cos 4 = -0.6536436$ , $\sin 4 = -0.7568025$ , $e^{-2} = 0.1353353$ , $\sin 2 = 0.9092974$ and $e^{-1} = 0.3678794$ ,

g(2) = 0.2279344 , \qquad 3k(1) = 0.5017677 , \qquad y(3) = 0.7297022 .

Differentiating gives $y'(t) = H(t-1)k(t-1) + 3H(t-2)k'(t-2)$ with $k'(t) = \tfrac12 e^{-t}(2\cos 2t - \sin 2t)$ , so

y'(3) = k(2) + 3k'(1) = -0.0512110 - 0.9610433 = -1.0122544 :

three units of time in, the blow has pushed the oscillator past its static level and it is on its way back down. Integrating the equation numerically from $t = 0$ to $t = 3$ , with the impulse imposed as a jump of $3$ in $y'$ at $t = 2$ , reproduces both numbers to seven decimals.

6. What the transform saved. Done by undetermined coefficients, this problem is three separate problems on $[0,1]$ , $[1,2]$ and $[2,\infty)$ with two matching conditions at each junction, one of which is a jump condition. The transform never splits the time axis: the three switch-on times appear as the three factors $1, e^{-s}, e^{-2s}$ , and the shifting theorem puts them back at the end.

Example(A repeated pole is resonance)Original

Solve $y'' + \omega^2 y = \sin \omega t$ with $y(0) = y'(0) = 0$ and $\omega > 0$ : an undamped oscillator driven exactly at its own frequency.

1. Transform. $p(s) = s^2 + \omega^2$ and $\mathcal{L}\{\sin\omega t\} = \omega/(s^2+\omega^2)$ , so

Y(s) = \frac{\omega}{(s^2+\omega^2)^2} .

The forcing has contributed poles at $s = \pm i\omega$ that the equation already had, so the poles are double. Partial fractions cannot reduce a double pole to simple ones, and that is the algebraic signature of resonance.

2. Invert using transform differentiation (step 5 of the method). From $\mathcal{L}\{\sin\omega t\} = \omega/(s^2+\omega^2)$ and $\mathcal{L}\{\cos\omega t\} = s/(s^2+\omega^2)$ , the rule $\mathcal{L}\{t\,v(t)\}(s) = -V'(s)$ gives

\mathcal{L}\{t\sin\omega t\} = \frac{2\omega s}{(s^2+\omega^2)^2} , \qquad \mathcal{L}\{t\cos\omega t\} = -\frac{d}{ds}\frac{s}{s^2+\omega^2} = \frac{s^2 - \omega^2}{(s^2+\omega^2)^2} .

Subtract the second from $\mathcal{L}\{\sin\omega t\}$ after matching coefficients:

\frac{\omega}{s^2+\omega^2} - \omega\,\frac{s^2-\omega^2}{(s^2+\omega^2)^2} = \frac{\omega(s^2+\omega^2) - \omega s^2 + \omega^3}{(s^2+\omega^2)^2} = \frac{2\omega^3}{(s^2+\omega^2)^2} ,

so $Y(s) = \dfrac{\omega}{(s^2+\omega^2)^2} = \dfrac{1}{2\omega^2}\left[\dfrac{\omega}{s^2+\omega^2} - \omega\,\dfrac{s^2-\omega^2}{(s^2+\omega^2)^2}\right]$ and

y(t) = \frac{1}{2\omega^2}\bigl(\sin\omega t - \omega t\cos\omega t\bigr) .

3. Checks. $y(0) = 0$ ; differentiating, $y' = \frac{1}{2\omega^2}(\omega\cos\omega t - \omega\cos\omega t + \omega^2 t \sin\omega t) = \tfrac12 t\sin\omega t$ , so $y'(0) = 0$ ; differentiating again, $y'' = \tfrac12\sin\omega t + \tfrac12 \omega t\cos\omega t$ , and

y'' + \omega^2 y = \tfrac12 \sin\omega t + \tfrac12\omega t \cos\omega t + \tfrac12 \sin\omega t - \tfrac12 \omega t\cos\omega t = \sin\omega t .

4. A number. At $\omega = 3$ , $t = \pi$ : $\sin 3\pi = 0$ and $\cos 3\pi = -1$ , so $y(\pi) = \frac{1}{18}(0 + 3\pi) = \pi/6 = 0.5235988$ .

5. Reading it. The amplitude $\frac{t}{2\omega}$ of the $\cos\omega t$ term grows without bound: an undamped oscillator driven at resonance absorbs energy indefinitely. Compare the damped resonance of nonhomogeneous systems, where the same double pole occurs but sits strictly left of the imaginary axis, so the factor $t$ is beaten by the decaying exponential and the response stays bounded. The dividing line is exactly whether the repeated root of $p$ has negative real part, which is why a vanishing damping coefficient is a singular limit and not a small perturbation.

ApplicationOriginal

Circuit transients. A series $RLC$ loop driven by a switch obeys $L\ddot q + R \dot q + q/C = V(t)$ , and every practical $V$ — a switch closing, a square wave, a spark — is a combination of steps and impulses. The transform turns the whole analysis into algebra in $s$ : the impedance $Ls + R + 1/(Cs)$ is $p(s)$ up to a factor, the transfer function $1/p(s)$ is what an engineer reads off a Bode plot, and the poles of $1/p$ are the damped natural frequencies. That the same formalism describes the mechanical oscillator is the dictionary $L \leftrightarrow m$ , $R \leftrightarrow b$ , $1/C \leftrightarrow k$ that Arfken sets out as the $RLC$ analogue ((Arfken, Weber & Harris (2013), §20.8, Eqs. (20.162)-(20.163), p. 1022)).
Impulse response and transfer function as the language of linear systems. Because $Y = F/p$ for zero initial data, a linear time-invariant system is completely described either by $k(t)$ in the time domain or by $1/p(s)$ in the transform domain, and measuring one determines the other. Stability is then a statement about where the poles are — all in the open left half plane — which is why control engineering is done with pole diagrams and why the transform variable $s$ carries the name complex frequency.
Instrument calibration from a kick. A ballistic galvanometer, a seismometer and an accelerometer are all read by delivering a brief pulse and recording the response. The statement above says why the recorded curve is the impulse response shifted to the arrival time, and why the area of the pulse — not its shape — sets the amplitude: the only thing a short pulse does is jump the highest derivative by its integral, exactly as Arfken's impulse transfers $P$ units of momentum ((Arfken, Weber & Harris (2013), §20.8, Example 20.8.4, p. 1019)).
Boundary data for partial differential equations. Transforming in $t$ reduces an evolution equation in one space variable to an ordinary differential equation in that variable, with the time-dependent boundary data appearing as transformed boundary values; inverting at the end often needs nothing but the shifting theorem. That is how a signal injected at the end of a transmission line comes out as $H(t - x/v)E(0, t-x/v)$ , a wave travelling at speed $v$ with nothing ahead of its front ((Arfken, Weber & Harris (2013), §20.8, Example 20.8.6, pp. 1023-1024); (Strauss (2008), §12.5, Examples 4 and 5, pp. 355-356)), and the same route handles the half-line heat equation with a prescribed surface temperature.
Why numerical work still needs the pole picture. A stiff constant-coefficient problem is one whose $p$ has roots of very different magnitude, and the stability of an explicit integrator is decided by where the products (step size) $\times$ (root) lie. The transform is thus the analytic companion to stiffness: it says what the exact solution does, which is what a numerical method has to be asked to reproduce.