#set document(title: "3.8 Matrix exponentials", author: "Jiří Lebl")
#set page(width: 8.5in, height: auto, margin: 1in)
#import "@preview/cetz:0.5.2"
#set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en")
#show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math"))
#set par(justify: true, leading: 0.62em, spacing: 0.9em)
#set enum(spacing: 1.1em)   // room between list items so tall inline fractions don't collide
#set list(spacing: 1.1em)
#set table(stroke: 0.5pt + rgb("#c7ccd3"))

#let BLUE = rgb("#183B6F")     // brand navy — section bars + example/solution labels (white on navy 11.09:1)
#let ORANGE = rgb("#A94509")   // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text)
#let RED = rgb("#DC2626")      // brand error-600
#let GREEN = rgb("#059669")    // brand success-600 (decoration only; small green text uses green-text #007942)

#show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt,
  fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt,
  text(fill: white, weight: "bold", size: 19pt, it.body))
#show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt,
  fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt,
  text(fill: white, weight: "bold", size: 12pt, it.body))
#show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body)
#show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body)

#let examplebox(label, title, body) = block(width: 100%, breakable: true,
  fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt,
  above: 12pt, below: 12pt)[
  #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt,
    text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]]
  #body]

// rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade
#let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint,
    stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt),
    above: 11pt, below: 11pt)[
    #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)]
    #linebreak() #body]

#let solutionbox(body) = block(above: 4pt, below: 8pt)[
  #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body]

#let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"),
  stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[
  #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))]

// Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a
// consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the
// print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a
// degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero.
#let _STD_W = 3.5
#let _WIDE_W = 5.6
#let _MAX_H = 3.4
#let _ASPECT_WIDE = 2.2
#let _UPSCALE_MAX = 1.15
#let stdfig(body) = context {
  let m = measure(body)
  let w = calc.max(m.width / 1in, 0.01)
  let h = calc.max(m.height / 1in, 0.01)
  let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W }
  let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX)
  align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body)))
}

#show figure: set block(breakable: false)
#set figure(gap: 8pt)
#show figure.caption: set text(size: 8.5pt, fill: rgb("#555"))


    
== 3.8#h(0.6em)Matrix exponentials


    
=== Definition


    In this section we present a different way of finding the fundamental matrix solution of a system. Suppose that we have the constant coefficient equation


    #math.equation(block: true, alt: "x → prime equals P x →")[$arrow(x)^(′) = P arrow(x)$]


    as usual. Now suppose that this was one equation (#math.equation(block: false, alt: "P")[$P$] is a number or a #math.equation(block: false, alt: "1 times 1")[$1 × 1$] matrix). Then the solution to this would be


    #math.equation(block: true, alt: "x → equals e to the power P t .")[$arrow(x) = e^(P t) .$]


    That doesn’t make sense if #math.equation(block: false, alt: "P")[$P$] is a larger matrix, but essentially the same computation that led to the above works for matrices when we define #math.equation(block: false, alt: "e to the power P t")[$e^(P t)$] properly. First let us write down the Taylor series for #math.equation(block: false, alt: "e to the power a t")[$e^(a t)$]for some number #math.equation(block: false, alt: "a")[$a$].


    #math.equation(block: true, alt: "e to the power a t equals 1 plus a t plus the fraction open parenthesis a t close parenthesis squared over 2 plus the fraction open parenthesis a t close parenthesis cubed over 6 plus the fraction open parenthesis a t close parenthesis to the power 4 over 24 plus ⋯ equals ∑ k equals 0 infinity the fraction open parenthesis a t close parenthesis to the power k over k !")[$e^(a t) = 1 + a t + frac(( a t )^(2), 2) + frac(( a t )^(3), 6) + frac(( a t )^(4), 24) + ⋯ = ∑_(k = 0)^(∞) frac(( a t )^(k), k !)$]


    Recall #math.equation(block: false, alt: "k ! equals 1 times 2 times 3 ⋯ k")[$k ! = 1 · 2 · 3 ⋯ k$] is the factorial, and #math.equation(block: false, alt: "0 ! equals 1")[$0 ! = 1$]. We differentiate this series term by term


    #math.equation(block: true, alt: "the fraction d over d t open parenthesis e to the power a t close parenthesis equals a plus a squared t plus the fraction a cubed t squared over 2 plus the fraction a to the power 4 t cubed over 6 plus ⋯ equals a open parenthesis 1 plus a t the fraction open parenthesis a t close parenthesis squared over 2 plus the fraction open parenthesis a t close parenthesis cubed over 6 plus ⋯ close parenthesis equals a e to the power a t .")[$frac(d, d t) ( e^(a t) ) = a + a^(2) t + frac(a^(3) t^(2), 2) + frac(a^(4) t^(3), 6) + ⋯ = a ( 1 + a t frac(( a t )^(2), 2) + frac(( a t )^(3), 6) + ⋯ ) = a e^(a t) .$]


    Maybe we can try the same trick with matrices. Suppose that for an #math.equation(block: false, alt: "n times n")[$n × n$] matrix #math.equation(block: false, alt: "A")[$A$] we define the #emph[matrix exponential] as


    #math.equation(block: true, alt: "e to the power A equals def 𝐼 plus A plus the fraction 1 over 2 A squared plus the fraction 1 over 6 A cubed plus ⋯ plus the fraction 1 over k ! A to the power k plus ⋯")[$e^(A) limits(=)^("def") 𝐼 + A + frac(1, 2) A^(2) + frac(1, 6) A^(3) + ⋯ + frac(1, k !) A^(k) + ⋯$]


    Let us not worry about convergence. The series really does always converge. We usually write #math.equation(block: false, alt: "P t")[$P t$]as #math.equation(block: false, alt: "t P")[$t P$] by convention when #math.equation(block: false, alt: "P")[$P$] is a matrix. With this small change and by the exact same calculation as above we have that


    #math.equation(block: true, alt: "the fraction d over d t open parenthesis e to the power t P close parenthesis equals P e to the power t P .")[$frac(d, d t) ( e^(t P) ) = P e^(t P) .$]


    Now #math.equation(block: false, alt: "P")[$P$] and hence #math.equation(block: false, alt: "e to the power t P")[$e^(t P)$] is an #math.equation(block: false, alt: "n times n")[$n × n$] matrix. What we are looking for is a vector. We note that in the #math.equation(block: false, alt: "1 times 1")[$1 × 1$] case we would at this point multiply by an arbitrary constant to get the general solution. In the matrix case we multiply by a column vector #math.equation(block: false, alt: "c →")[$arrow(c)$].


    #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
Let #math.equation(block: false, alt: "P")[$P$] be an #math.equation(block: false, alt: "n times n")[$n × n$] matrix. Then the general solution to #math.equation(block: false, alt: "x → prime equals P x →")[$arrow(x)^(′) = P arrow(x)$] is


      #math.equation(block: true, alt: "x → equals e to the power t P c → ,")[$arrow(x) = e^(t P) arrow(c) ,$]


      where #math.equation(block: false, alt: "c →")[$arrow(c)$] is an arbitrary constant vector. In fact #math.equation(block: false, alt: "x → open parenthesis 0 close parenthesis equals c →")[$arrow(x) ( 0 ) = arrow(c)$].
]


    Let us check.


    #math.equation(block: true, alt: "the fraction d over d t x → equals the fraction d over d t open parenthesis e to the power t P c → close parenthesis equals P e to the power t P c → equals P x → .")[$frac(d, d t) arrow(x) = frac(d, d t) ( e^(t P) arrow(c) ) = P e^(t P) arrow(c) = P arrow(x) .$]


    Hence #math.equation(block: false, alt: "e to the power t P")[$e^(t P)$]is the #emph[fundamental matrix solution] of the homogeneous system. If we find a way to compute the matrix exponential, we will have another method of solving constant coefficient homogeneous systems. It also makes it easy to solve for initial conditions. To solve #math.equation(block: false, alt: "x → prime equals A x →")[$arrow(x)^(′) = A arrow(x)$], #math.equation(block: false, alt: "x → open parenthesis 0 close parenthesis equals b →")[$arrow(x) ( 0 ) = arrow(b)$] we take the solution


    #math.equation(block: true, alt: "x → equals e to the power t A b →")[$arrow(x) = e^(t A) arrow(b)$]


    This equation follows because #math.equation(block: false, alt: "e to the power 0 A equals 𝐼")[$e^(0 A) = 𝐼$], so #math.equation(block: false, alt: "x → open parenthesis 0 close parenthesis equals e to the power 0 A b → equals b →")[$arrow(x) ( 0 ) = e^(0 A) arrow(b) = arrow(b)$].


    We mention a drawback of matrix exponentials. In general #math.equation(block: false, alt: "e to the power A plus B not equal to e to the power A e to the power B")[$e^(A + B) ≠ e^(A) e^(B)$]. The trouble is that matrices do not commute, that is, in general #math.equation(block: false, alt: "A B not equal to B A")[$A B ≠ B A$]. If you try to prove #math.equation(block: false, alt: "e to the power A plus B not equal to e to the power A e to the power B")[$e^(A + B) ≠ e^(A) e^(B)$] using the Taylor series, you will see why the lack of commutativity becomes a problem. However, it is still true that if #math.equation(block: false, alt: "A B equals B A")[$A B = B A$], that is, if #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$] #emph[commute], then #math.equation(block: false, alt: "e to the power A plus B equals e to the power A e to the power B")[$e^(A + B) = e^(A) e^(B)$]. We will find this fact useful. Let us restate this as a theorem to make a point.


    #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
If #math.equation(block: false, alt: "A B equals B A")[$A B = B A$], then #math.equation(block: false, alt: "e to the power A plus B equals e to the power A e to the power B")[$e^(A + B) = e^(A) e^(B)$]. Otherwise #math.equation(block: false, alt: "e to the power A plus B not equal to e to the power A e to the power B")[$e^(A + B) ≠ e^(A) e^(B)$] in general.
]


    
=== Simple cases


    In some instances it may work to just plug into the series definition. Suppose the matrix is diagonal. For example, #math.equation(block: false, alt: "D equals [ a, 0; 0, b ]")[$D = [ a & 0 \ 0 & b ]$]. Then


    #math.equation(block: true, alt: "D to the power k equals [ a to the power k, 0; 0, b to the power k ]")[$D^(k) = [ a^(k) & 0 \ 0 & b^(k) ]$]


    and


    #math.equation(block: true, alt: "e to the power D equals 𝐼 plus D plus the fraction 1 over 2 D squared plus the fraction 1 over 6 D cubed plus ⋯; equals [ 1, 0; 0, 1 ] plus [ a, 0; 0, b ] plus the fraction 1 over 2 [ a squared, 0; 0, b squared ] plus the fraction 1 over 6 [ a cubed, 0; 0, b cubed ] plus ⋯ equals [ e to the power a, 0; 0, e to the power b ]")[$e^(D) = 𝐼 + D + frac(1, 2) D^(2) + frac(1, 6) D^(3) + ⋯ \ = [ 1 & 0 \ 0 & 1 ] + [ a & 0 \ 0 & b ] + frac(1, 2) [ a^(2) & 0 \ 0 & b^(2) ] + frac(1, 6) [ a^(3) & 0 \ 0 & b^(3) ] + ⋯ = [ e^(a) & 0 \ 0 & e^(b) ]$]


    So by this rationale we have that


    #math.equation(block: true, alt: "e to the power 𝐼 equals [ e, 0; 0, e ] and e to the power a 𝐼 equals [ e to the power a, 0; 0, e to the power a ]")[$e^(𝐼) = [ e & 0 \ 0 & e ] #h(1em) "and" #h(1em) e^(a 𝐼) = [ e^(a) & 0 \ 0 & e^(a) ]$]


    This makes exponentials of certain other matrices easy to compute. Notice for example that the matrix #math.equation(block: false, alt: "A equals [ 5, 4; minus 1, 1 ]")[$A = [ 5 & 4 \ − 1 & 1 ]$] can be written as #math.equation(block: false, alt: "3 I plus B")[$3 italic(I) + italic(B)$] where #math.equation(block: false, alt: "𝐵 equals [ 2, 4; minus 1, minus 2 ]")[$𝐵 = [ 2 & 4 \ − 1 & − 2 ]$]. Notice that #math.equation(block: false, alt: "B squared equals [ 0, 0; 0, 0 ]")[$italic(B)^(2) = [ 0 & 0 \ 0 & 0 ]$]. So #math.equation(block: false, alt: "B to the power k equals 0")[$italic(B)^(italic(k)) = 0$] for all #math.equation(block: false, alt: "k greater than or equal to 2")[$k ≥ 2$]. Therefore, #math.equation(block: false, alt: "e to the power B equals I plus B")[$e^(B) = italic(I) + italic(B)$]. Suppose we actually want to compute #math.equation(block: false, alt: "e to the power t A")[$e^(t A)$]. The matrices #math.equation(block: false, alt: "3 t I")[$3 italic(t) italic(I)$] and #math.equation(block: false, alt: "t B")[$t B$] commute (exercise: check this) and #math.equation(block: false, alt: "e to the power t B equals I plus t B")[$e^(t B) = italic(I) + italic(t) italic(B)$], since #math.equation(block: false, alt: "open parenthesis t B close parenthesis squared equals t squared B squared equals 0")[$attach(( t B ), t: 2) = t^(2) B^(2) = 0$]. We write


    #math.equation(block: true, alt: "e to the power t A equals e to the power 3 t I plus t 𝐵 equals e to the power 3 t 𝐼 e to the power t 𝐵 equals [ e to the power 3 t, 0; 0, e to the power 3 t ] open parenthesis I plus t B close parenthesis; equals [ e to the power 3 t, 0; 0, e to the power 3 t ] [ 1 plus 2 t, 4 t; minus t, 1 minus 2 t ] equals [ open parenthesis 1 plus 2 t close parenthesis e to the power 3 t, 4 t e to the power 3 t; minus t e to the power 3 t, open parenthesis 1 minus 2 t close parenthesis e to the power 3 t ]")[$e^(t A) = e^(3 italic(t) italic(I) + t 𝐵) = e^(3 t 𝐼) e^(t 𝐵) = [ e^(3 t) & 0 \ 0 & e^(3 t) ] ( I + t B ) \ = [ e^(3 t) & 0 \ 0 & e^(3 t) ] [ 1 + 2 t & 4 t \ − t & 1 − 2 t ] = [ ( 1 + 2 t ) e^(3 t) & 4 t e^(3 t) \ − t e^(3 t) & ( 1 − 2 t ) e^(3 t) ]$]


    So we have found the fundamental matrix solution for the system #math.equation(block: false, alt: "x → prime equals A x →")[$arrow(x)^(′) = A arrow(x)$]. Note that this matrix has a repeated eigenvalue with a defect; there is only one eigenvector for the eigenvalue 3. So we have found a perhaps easier way to handle this case. In fact, if a matrix #math.equation(block: false, alt: "A")[$A$] is #math.equation(block: false, alt: "2 times 2")[$2 × 2$] and has an eigenvalue #math.equation(block: false, alt: "λ")[$λ$] of multiplicity 2, then either #math.equation(block: false, alt: "A")[$A$] is diagonal, or #math.equation(block: false, alt: "A equals λ 𝐼 plus B")[$A = λ 𝐼 + B$] where #math.equation(block: false, alt: "B squared equals 0")[$B^(2) = 0$]. This is a good exercise.


    #notebox("Your Turn", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
Suppose that #math.equation(block: false, alt: "A")[$A$] is #math.equation(block: false, alt: "2 times 2")[$2 × 2$] and #math.equation(block: false, alt: "λ")[$λ$] is the only eigenvalue. Then show that #math.equation(block: false, alt: "open parenthesis A minus λ 𝐼 close parenthesis squared equals 0")[$attach(( A − λ 𝐼 ), t: 2) = 0$]. Then we can write #math.equation(block: false, alt: "A equals λ 𝐼 plus B")[$A = λ 𝐼 + B$], where #math.equation(block: false, alt: "B squared equals 0")[$B^(2) = 0$]. Hint: First write down what does it mean for the eigenvalue to be of multiplicity #math.equation(block: false, alt: "2")[$2$]. You will get an equation for the entries. Now compute the square of #math.equation(block: false, alt: "B")[$B$].
]


    Matrices #math.equation(block: false, alt: "B")[$B$] such that #math.equation(block: false, alt: "B to the power k equals 0")[$B^(k) = 0$] for some #math.equation(block: false, alt: "k")[$k$] are called #emph[nilpotent]. Computation of the matrix exponential for nilpotent matrices is easy by just writing down the first #math.equation(block: false, alt: "k")[$k$] terms of the Taylor series.


    
=== General Matrices


    In general, the exponential is not as easy to compute as above. We usually cannot write a matrix as a sum of commuting matrices where the exponential is simple for each one. But fear not, it is still not too difficult provided we can find enough eigenvectors. First we need the following interesting result about matrix exponentials. For two square matrices #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$], with #math.equation(block: false, alt: "B")[$B$] #emph[invertible], we have


    #math.equation(block: true, alt: "e to the power B A B to the power minus 1 equals B e to the power A B to the power minus 1 .")[$e^(B A B^(− 1)) = B e^(A) B^(− 1) .$]


    This can be seen by writing down the Taylor series. First note that


    #math.equation(block: true, alt: "open parenthesis B A B to the power minus 1 close parenthesis squared equals B A B to the power minus 1 B A B to the power minus 1 equals B A 𝐼 A B to the power minus 1 equals B A squared B to the power minus 1")[$attach(( B A B^(− 1) ), t: 2) = B A B^(− 1) B A B^(− 1) = B A 𝐼 A B^(− 1) = B A^(2) B^(− 1)$]


    And hence by the same reasoning #math.equation(block: false, alt: "open parenthesis B A B to the power minus 1 close parenthesis to the power k equals B A to the power k B to the power minus 1")[$attach(( B A B^(− 1) ), t: k) = B A^(k) B^(− 1)$]. Now write down the Taylor series for #math.equation(block: false, alt: "e to the power B A B to the power minus 1")[$e^(B A B^(− 1))$].


    #math.equation(block: true, alt: "e to the power B A B to the power minus 1 equals 𝐼 plus B A B to the power minus 1 plus the fraction 1 over 2 open parenthesis B A B to the power minus 1 close parenthesis squared plus the fraction 1 over 6 open parenthesis B A B to the power minus 1 close parenthesis cubed plus ⋯; equals B B to the power minus 1 plus B A B to the power minus 1 plus the fraction 1 over 2 B A squared B to the power minus 1 plus the fraction 1 over 6 B A cubed B to the power minus 1 plus ⋯; equals B open parenthesis 𝐼 plus A plus the fraction 1 over 2 A squared plus the fraction 1 over 6 A cubed plus ⋯ close parenthesis B to the power minus 1; equals B e to the power A B to the power minus 1 .")[$e^(B A B^(− 1)) = 𝐼 + B A B^(− 1) + frac(1, 2) attach(( B A B^(− 1) ), t: 2) + frac(1, 6) attach(( B A B^(− 1) ), t: 3) + ⋯ \ = B B^(− 1) + B A B^(− 1) + frac(1, 2) B A^(2) B^(− 1) + frac(1, 6) B A^(3) B^(− 1) + ⋯ \ = B ( 𝐼 + A + frac(1, 2) A^(2) + frac(1, 6) A^(3) + ⋯ ) B^(− 1) \ = B e^(A) B^(− 1) .$]


    Given a square matrix #math.equation(block: false, alt: "A")[$A$], we can sometimes write #math.equation(block: false, alt: "A equals E D E to the power minus 1")[$A = E D E^(− 1)$], where #math.equation(block: false, alt: "D")[$D$] is diagonal and #math.equation(block: false, alt: "E")[$E$] invertible. This procedure is called #emph[diagonalization]. If we can do that, the computation of the exponential becomes easy. Adding #math.equation(block: false, alt: "t")[$t$] into the mix we see that we can then easily compute the exponential


    #math.equation(block: true, alt: "e to the power t A equals E e to the power t D E to the power minus 1 .")[$e^(t A) = E e^(t D) E^(− 1) .$]


    To diagonalize #math.equation(block: false, alt: "A")[$A$] we will need #math.equation(block: false, alt: "n")[$n$] #emph[linearly independent]eigenvectors of #math.equation(block: false, alt: "A")[$A$]. Otherwise this method of computing the exponential does not work and we need to be trickier, but we will not get into such details. We let #math.equation(block: false, alt: "E")[$E$] be the matrix with the eigenvectors as columns. Let #math.equation(block: false, alt: "λ sub 1 , λ sub 2 , ⋯ , λ sub n")[$λ_(1) , #h(0.222em) λ_(2) , ⋯ , λ_(n)$] be the eigenvalues and let #math.equation(block: false, alt: "v → sub 1 , v sub 2 → , ⋯ , v → sub n")[$arrow(v)_(1) , #h(0.222em) arrow(v_(2)) , ⋯ , arrow(v)_(n)$] be the eigenvectors, then #math.equation(block: false, alt: "E equals [ v → sub 1 v → sub 2 ⋯ v → sub n ]")[$E = [ arrow(v)_(1) " " " " arrow(v)_(2) " " " " ⋯ " " " " arrow(v)_(n) ]$]. Let #math.equation(block: false, alt: "D")[$D$] be the diagonal matrix with the eigenvalues on the main diagonal. That is


    #math.equation(block: true, alt: "D equals [ λ sub 1, 0, ⋯, 0; 0, λ sub 2, ⋯, 0; ⋮, ⋮, ⋱, ⋮; 0, 0, ⋯, λ sub n ]")[$D = [ λ_(1) & 0 & ⋯ & 0 \ 0 & λ_(2) & ⋯ & 0 \ ⋮ & ⋮ & ⋱ & ⋮ \ 0 & 0 & ⋯ & λ_(n) ]$]


    We compute


    #math.equation(block: true, alt: "A E equals A [ v → sub 1, v → sub 2, ⋯, v → sub n ]; equals [ A v → sub 1, A v → sub 2, ⋯, A v → sub 3 ]; equals [ λ sub 1 v → sub 1, λ sub 2 v → sub 2, ⋯, λ sub n v → sub n ]; equals [ v → sub 1, v → sub 2, ⋯, v → sub n ] D; equals E D .")[$A E = A [ arrow(v)_(1) & arrow(v)_(2) & ⋯ & arrow(v)_(n) ] \ = [ A arrow(v)_(1) & A arrow(v)_(2) & ⋯ & A arrow(v)_(3) ] \ = [ λ_(1) arrow(v)_(1) & λ_(2) arrow(v)_(2) & ⋯ & λ_(n) arrow(v)_(n) ] \ = [ arrow(v)_(1) & arrow(v)_(2) & ⋯ & arrow(v)_(n) ] D \ = E D .$]


    The columns of #math.equation(block: false, alt: "E")[$E$] are linearly independent as these are linearly independent eigenvectors of #math.equation(block: false, alt: "A")[$A$]. Hence #math.equation(block: false, alt: "E")[$E$] is #emph[invertible]. Since #math.equation(block: false, alt: "A E equals E D")[$A E = E D$], we right multiply by #math.equation(block: false, alt: "E to the power minus 1")[$E^(− 1)$] and we get


    #math.equation(block: true, alt: "A equals E D E to the power minus 1 .")[$A = E D E^(− 1) .$]


    This means that . #math.equation(block: false, alt: "e to the power A equals E e to the power D E to the power minus 1")[$e^(A) = E e^(D) E^(− 1)$] Multiplying the matrix by #math.equation(block: false, alt: "t")[$t$] we obtain


    #math.equation(block: true, alt: "e to the power t A equals E e to the power t D E to the power minus 1 equals E [ e to the power λ sub 1 t, 0, ⋯, 0; 0, e to the power λ sub 2 t, ⋯, 0; ⋮, ⋮, ⋱, ⋮; 0, 0, ⋯, e to the power λ sub n t ] E to the power minus 1")[$e^(t A) = E e^(t D) E^(− 1) = E [ e^(λ_(1) t) & 0 & ⋯ & 0 \ 0 & e^(λ_(2) t) & ⋯ & 0 \ ⋮ & ⋮ & ⋱ & ⋮ \ 0 & 0 & ⋯ & e^(λ_(n) t) ] E^(− 1)$]


    The formula (3.8.1), therefore, gives the formula for computing the fundamental matrix solution #math.equation(block: false, alt: "e to the power t A")[$e^(t A)$] for the system #math.equation(block: false, alt: "x → prime equals A x →")[$arrow(x)^(′) = A arrow(x)$], in the case where we have #math.equation(block: false, alt: "n")[$n$] linearly independent eigenvectors.


    Notice that this computation still works when the eigenvalues and eigenvectors are complex, though then you will have to compute with complex numbers. It is clear from the definition that if #math.equation(block: false, alt: "A")[$A$] is real, then #math.equation(block: false, alt: "e to the power t A")[$e^(t A)$] is real. So you will only need complex numbers in the computation and you may need to apply #link("http://www.jirka.org/diffyqs/htmlver/diffyqsse11.html#x16-36005r2")[Euler’s formula] to simplify the result. If simplified properly the final matrix will not have any complex numbers in it.


    #examplebox("Example 1")[][
Compute the fundamental matrix solution using the matrix exponentials for the system


        #math.equation(block: true, alt: "[ to the power x; y equals [ 1, 2; 2, 1 ] [ x; y ] .")[$\[^(x \ y) = \[ 1 & 2 \ 2 & 1 \] \[ x \ y \] .$]


        Then compute the particular solution for the initial conditions #math.equation(block: false, alt: "x open parenthesis 0 close parenthesis equals 4")[$x ( 0 ) = 4$] and #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 2")[$y ( 0 ) = 2$].


        Let #math.equation(block: false, alt: "A")[$A$] be the coefficient matrix #math.equation(block: false, alt: "[ 1, 2; 2, 1 ]")[$[ 1 & 2 \ 2 & 1 ]$]. We first compute (exercise) that the eigenvalues are 3 and -1 and corresponding eigenvectors are #math.equation(block: false, alt: "[ 1; 1 ]")[$[ 1 \ 1 ]$] and #math.equation(block: false, alt: "[ 1; minus 1 ]")[$[ 1 \ − 1 ]$]. Hence the diagonalization of #math.equation(block: false, alt: "A")[$A$] is


        #math.equation(block: true, alt: "[ 1, 2; 2, 1 ] ⏟ A equals [ 1, 1; 1, minus 1 ] ⏟ E [ 3, 0; 0, minus 1 ] ⏟ D [ to the power 1, 1; 1, minus 1 ⏟ E to the power minus 1 .")[$limits(underbrace(\[ 1 & 2 \ 2 & 1 \]))_(A) = limits(underbrace(\[ 1 & 1 \ 1 & − 1 \]))_(E) limits(underbrace(\[ 3 & 0 \ 0 & − 1 \]))_(D) limits(underbrace(\[^(1 & 1 \ 1 & − 1)))_(E^(− 1)) .$]


        We write


        #math.equation(block: true, alt: "e to the power t A equals E e to the power t D E to the power minus 1 equals [ 1, 1; 1, minus 1 ] [ e to the power 3 t, 0; 0, e to the power minus t ] [ to the power 1, 1; 1, minus 1; equals [ 1, 1; 1, minus 1 ] [ e to the power 3 t, 0; 0, e to the power minus t ] the fraction minus 1 over 2 [ minus 1, minus 1; minus 1, 1 ]; equals the fraction minus 1 over 2 [ e to the power 3 t, e to the power minus t; e to the power 3 t, minus e to the power minus t ] [ minus 1, minus 1; minus 1, 1 ]; equals the fraction minus 1 over 2 [ minus e to the power 3 t minus e to the power minus t, minus e to the power 3 t plus e to the power minus t; minus e to the power 3 t plus e to the power minus t, minus e to the power 3 t minus e to the power minus t ] equals [ the fraction e to the power 3 t plus e to")[$e^(t A) = E e^(t D) E^(− 1) = \[ 1 & 1 \ 1 & − 1 \] \[ e^(3 t) & 0 \ 0 & e^(− t) \] \[^(1 & 1 \ 1 & − 1) \ = \[ 1 & 1 \ 1 & − 1 \] \[ e^(3 t) & 0 \ 0 & e^(− t) \] frac(− 1, 2) \[ − 1 & − 1 \ − 1 & 1 \] \ = frac(− 1, 2) \[ e^(3 t) & e^(− t) \ e^(3 t) & − e^(− t) \] \[ − 1 & − 1 \ − 1 & 1 \] \ = frac(− 1, 2) \[ − e^(3 t) − e^(− t) & − e^(3 t) + e^(− t) \ − e^(3 t) + e^(− t) & − e^(3 t) − e^(− t) \] = \[ frac(e^(3 t) + e^(− t), 2) & frac(e^(3 t) − e^(− t), 2) \ frac(e^(3 t) − e^(− t), 2) & frac(e^(3 t) + e^(− t), 2) \] .$]


        The initial conditions are #math.equation(block: false, alt: "x open parenthesis 0 close parenthesis equals 4")[$x ( 0 ) = 4$] and #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 2")[$y ( 0 ) = 2$]. Hence, by the property that #math.equation(block: false, alt: "e to the power 0 A equals 𝐼")[$e^(0 A) = 𝐼$] we find that the particular solution we are looking for is #math.equation(block: false, alt: "e to the power t A b →")[$e^(t A) arrow(b)$] where #math.equation(block: false, alt: "b →")[$arrow(b)$] is #math.equation(block: false, alt: "[ 4; 2 ]")[$[ 4 \ 2 ]$]. Then the particular solution we are looking for is


        #math.equation(block: true, alt: "[ x; y ] equals [ the fraction e to the power 3 t plus e to the power minus t over 2, the fraction e to the power 3 t minus e to the power minus t over 2; the fraction e to the power 3 t minus e to the power minus t over 2, the fraction e to the power 3 t plus e to the power minus t over 2 ] [ 4; 2 ] equals [ 2 e to the power 3 t plus 2 e to the power minus t plus e to the power 3 t minus e to the power minus t; 2 e to the power 3 t minus 2 e to the power minus t plus e to the power 3 t plus e to the power minus t ] equals [ 3 e to the power 3 t plus e to the power minus t; 3 e to the power 3 t minus e to the power minus t ]")[$[ x \ y ] = [ frac(e^(3 t) + e^(− t), 2) & frac(e^(3 t) − e^(− t), 2) \ frac(e^(3 t) − e^(− t), 2) & frac(e^(3 t) + e^(− t), 2) ] [ 4 \ 2 ] = [ 2 e^(3 t) + 2 e^(− t) + e^(3 t) − e^(− t) \ 2 e^(3 t) − 2 e^(− t) + e^(3 t) + e^(− t) ] = [ 3 e^(3 t) + e^(− t) \ 3 e^(3 t) − e^(− t) ]$]
]


    
=== Fundamental Matrix Solutions


    We note that if you can compute the fundamental matrix solution in a different way, you can use this to find the matrix exponential #math.equation(block: false, alt: "e to the power t A")[$e^(t A)$]. The fundamental matrix solution of a system of ODEs is not unique. The exponential is the fundamental matrix solution with the property that for #math.equation(block: false, alt: "t equals 0")[$t = 0$] we get the identity matrix. So we must find the right fundamental matrix solution. Let #math.equation(block: false, alt: "X")[$X$] be any fundamental matrix solution to #math.equation(block: false, alt: "x → prime equals A x →")[$arrow(x)^(′) = A arrow(x)$]. Then we claim


    #math.equation(block: true, alt: "e to the power t A equals X open parenthesis t close parenthesis [ X open parenthesis 0 close parenthesis ] to the power minus 1 .")[$e^(t A) = X ( t ) [ X ( 0 ) ]^(− 1) .$]


    Clearly, if we plug #math.equation(block: false, alt: "t equals 0")[$t = 0$] into #math.equation(block: false, alt: "X open parenthesis t close parenthesis [ X open parenthesis 0 close parenthesis ] to the power minus 1")[$X ( t ) [ X ( 0 ) ]^(− 1)$] we get the identity. We can multiply a fundamental matrix solution on the right by any constant invertible matrix and we still get a fundamental matrix solution. All we are doing is changing what the arbitrary constants are in the general solution #math.equation(block: false, alt: "x → open parenthesis t close parenthesis equals X open parenthesis t close parenthesis c →")[$arrow(x) ( t ) = X ( t ) arrow(c)$].


    
=== Approximations


    If you think about it, the computation of any fundamental matrix solution #math.equation(block: false, alt: "X")[$X$] using the eigenvalue method is just as difficult as the computation of #math.equation(block: false, alt: "e to the power t A")[$e^(t A)$]. So perhaps we did not gain much by this new tool. However, the Taylor series expansion actually gives us a very easy way to approximate solutions, which the eigenvalue method did not.


    The simplest thing we can do is to just compute the series up to a certain number of terms. There are better ways to approximate the exponential#math.equation(block: false, alt: "to the power 1")[$1$]. In many cases however, few terms of the Taylor series give a reasonable approximation for the exponential and may suffice for the application. For example, let us compute the first #math.equation(block: false, alt: "4")[$4$] terms of the series for the matrix #math.equation(block: false, alt: "1 2; 2 1")[$1 2 \ 2 1$].


    #math.equation(block: true, alt: "e to the power t A approximately equals I plus t A plus the fraction t squared over 2 A squared plus the fraction t cubed over 6 A cubed equals I plus t [ 1, 2; 2, 1 ] plus t squared [ the fraction 5 over 2, 2; 2, the fraction 5 over 2 ] plus t cubed [ the fraction 13 over 6, the fraction 7 over 3; the fraction 7 over 3, the fraction 13 over 6 ] equals; equals [ 1 plus t plus the fraction 5 over 2 t squared plus the fraction 13 over 6 t cubed, 2 t plus 2 t squared plus the fraction 7 over 3 t cubed; 2 t plus 2 t squared plus the fraction 7 over 3 t cubed, 1 plus t plus the fraction 5 over 2 t squared plus the fraction 13 over 6 t cubed ] .")[$e^(t A) ≈ I + t A + frac(t^(2), 2) A^(2) + frac(t^(3), 6) A^(3) = I + t [ 1 & 2 \ 2 & 1 ] + t^(2) [ frac(5, 2) & 2 \ 2 & frac(5, 2) ] + t^(3) [ frac(13, 6) & frac(7, 3) \ frac(7, 3) & frac(13, 6) ] = \ = [ 1 + t + frac(5, 2) t^(2) + frac(13, 6) t^(3) & 2 t + 2 t^(2) + frac(7, 3) t^(3) \ 2 t + 2 t^(2) + frac(7, 3) t^(3) & 1 + t + frac(5, 2) t^(2) + frac(13, 6) t^(3) ] .$]


    Just like the scalar version of the Taylor series approximation, the approximation will be better for small #math.equation(block: false, alt: "t")[$t$] and worse for larger #math.equation(block: false, alt: "t")[$t$]. For larger #math.equation(block: false, alt: "t")[$t$], we will generally have to compute more terms. Let us see how we stack up against the real solution with #math.equation(block: false, alt: "t equals 0.1")[$t = 0.1$]. The approximate solution is approximately (rounded to #math.equation(block: false, alt: "8")[$8$] decimal places)


    #math.equation(block: true, alt: "e to the power 0.1 A approximately equals 𝐼 plus 0.1 A plus the fraction 0.1 squared over 2 plus the fraction 0.1 cubed over 6 A cubed equals [ 1.12716667, 0.22233333; 0.22233333, 1.12716667 ]")[$e^(0.1 A) ≈ 𝐼 + 0.1 A + frac(0.1^(2), 2) + frac(0.1^(3), 6) A^(3) = [ 1.12716667 & 0.22233333 \ 0.22233333 & 1.12716667 ]$]


    And plugging #math.equation(block: false, alt: "t equals 0.1")[$t = 0.1$] into the real solution (rounded to #math.equation(block: false, alt: "8")[$8$] decimal places) we get


    #math.equation(block: true, alt: "e to the power 0.1 A equals [ 1.12734811, 0.22251069; 0.22251069, 1.12734811 ]")[$e^(0.1 A) = [ 1.12734811 & 0.22251069 \ 0.22251069 & 1.12734811 ]$]


    Not bad at all! Although if we take the same approximation for #math.equation(block: false, alt: "t equals 1")[$t = 1$] we get


    #math.equation(block: true, alt: "𝐼 plus A plus the fraction 1 over 2 A squared plus the fraction 1 over 6 A cubed equals [ 6.66666667, 6.33333333; 6.33333333, 6.66666667 ]")[$𝐼 + A + frac(1, 2) A^(2) + frac(1, 6) A^(3) = [ 6.66666667 & 6.33333333 \ 6.33333333 & 6.66666667 ]$]


    while the real value is (again rounded to #math.equation(block: false, alt: "8")[$8$] decimal places)


    #math.equation(block: true, alt: "e to the power A equals [ 10.22670818, 9.85882874; 9.85882874, 10.22670818 ]")[$e^(A) = [ 10.22670818 & 9.85882874 \ 9.85882874 & 10.22670818 ]$]


    So the approximation is not very good once we get up to #math.equation(block: false, alt: "t equals 1")[$t = 1$]. To get a good approximation at #math.equation(block: false, alt: "t equals 1")[$t = 1$] (say up to #math.equation(block: false, alt: "2")[$2$] decimal places) we would need to go up to the #math.equation(block: false, alt: "11 to the power t h")[$11^(t h)$] power (exercise).


    
=== Footnotes


    \[1\] C. Moler and C.F. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later, SIAM Review 45 (1), 2003, 3–49