Chapter 8 Experiment 7: Eigenvalue Multiplicity and Dynamical Systems

This experiment is organized around the theme “one eigenvalue, two fates.” Through difference equations, differential equations, and the matrix exponential, it brings the multiplicity theory of §8.2 and the similarity transformations of §8.3–§8.4 to life in settings that can be computed and visualized.

Design of the Experiment: Two Contrasting Matrices¶

The whole experiment revolves around the following two $2 \times 2$ matrices:

\mathbf{A}_{\text{I}} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix}, \qquad \mathbf{A}_{\text{II}} = \begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix}

(1)

# ============================================================
# Environment setup: import the libraries and define the two featured matrices
# ============================================================
import numpy as np
import scipy.linalg                        # provides expm (the matrix exponential)
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
from matplotlib.gridspec import GridSpec

# --- Font settings ---
# Fonts: the English edition needs no CJK font, so nothing is downloaded when _lang == 'en';
# the Chinese editions use this same block to fetch Noto Sans TC/SC where no CJK font is installed
import os, urllib.request
import matplotlib.font_manager as fm
_lang = 'en'
_cjk = ['Microsoft JhengHei', 'PingFang TC', 'Noto Sans CJK TC', 'Noto Sans TC']
_have = {f.name for f in fm.fontManager.ttflist}
if _lang != 'en' and not _have & set(_cjk):
    _font = os.path.join(os.path.expanduser('~'), '.cache', 'fonts', 'NotoSansTC.ttf')
    try:
        if not os.path.exists(_font):
            os.makedirs(os.path.dirname(_font), exist_ok=True)
            urllib.request.urlretrieve('https://github.com/google/fonts/raw/main/ofl/notosanstc/NotoSansTC%5Bwght%5D.ttf', _font + '.part')
            os.replace(_font + '.part', _font)
        fm.fontManager.addfont(_font)
        _have.add('Noto Sans TC')
    except OSError as err:
        print('Could not download the CJK font; Chinese text in figures may not display:', err)
plt.rcParams['font.family'] = [f for f in _cjk if f in _have] + ['DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# --- book colors ---
C_VIOLET  = '#57068C'   # violet: Matrix I (diagonalizable)
C_BLUE    = '#006385'   # blue: Matrix II (not diagonalizable)
C_ORANGE  = '#FF5D47'   # orange: highlighting the contrast
C_TEAL    = '#2AD2C9'   # teal: auxiliary
LIGHT_BG    = '#F9F9FB'

# --- Display helper functions ---
def print_header(title):
    bar = '=' * (len(title) + 4)
    print(f'\n{bar}')
    print(f'  {title}')
    print(f'{bar}')

def print_step(n, desc):
    print(f'\n  Step {n}: {desc}')
    print('  ' + '-' * (len(desc) + 4))

def compare_print(label, val_a, val_b):
    print(f'  {label}')
    print(f'    Matrix I (diagonalizable): {val_a}')
    print(f'    Matrix II (Jordan block):  {val_b}')

# --- Define the two featured matrices ---
A_diag   = np.array([[2., 0.],             # Matrix I: purely diagonal, g=a=2
                     [0., 2.]])
A_jordan = np.array([[2., 1.],             # Matrix II: Jordan block, g=1<a=2
                     [0., 2.]])

print_header('The two featured matrices')
compare_print('The matrices themselves:',
    f'\n    {A_diag[0]}\n    {A_diag[1]}',
    f'\n    {A_jordan[0]}\n    {A_jordan[1]}')

# --- Check the geometric multiplicities ---
print_step('0', 'Check the geometric multiplicity (dimension of the eigenspace)')
for name, A in [('I', A_diag), ('II', A_jordan)]:
    null_mat = A - 2 * np.eye(2)           # (A - 2I)
    # use the SVD to compute the dimension of the null space (= 2 - rank)
    rank = np.linalg.matrix_rank(null_mat)
    geo_mult = 2 - rank
    print(f'  Matrix {name}: rank(A-2I) = {rank}, geometric multiplicity g_2 = {geo_mult}')

Part 1: The Difference Equation $\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k$ ¶

1.1 Computing $\mathbf{A}^k$ by Hand: The Key Role of the Binomial Expansion¶

Both matrices can be written in the form $\mathbf{A} = 2\mathbf{I} + \mathbf{N}$ ; the difference lies in $\mathbf{N}$ :

\mathbf{A}_{\text{I}} = 2\mathbf{I} + \mathbf{0}, \qquad \mathbf{A}_{\text{II}} = 2\mathbf{I} + \underbrace{\begin{pmatrix}0 & 1 \\ 0 & 0\end{pmatrix}}_{\mathbf{N},\ \mathbf{N}^2=\mathbf{0}}

(2)

Since $2\mathbf{I}$ commutes with every matrix, we may use the binomial expansion:

\mathbf{A}^k = (2\mathbf{I} + \mathbf{N})^k = \sum_{j=0}^{k} \binom{k}{j} (2\mathbf{I})^{k-j} \mathbf{N}^j

(3)

Why Does the Expansion Terminate at

j=1

?

For Matrix I, $\mathbf{N} = \mathbf{0}$ , so every term with $j \geq 1$ is zero and only $2^k \mathbf{I}$ remains.

For Matrix II, $\mathbf{N}^2 = \mathbf{0}$ (check: $\begin{pmatrix}0&1\\0&0\end{pmatrix}^2 = \mathbf{0}$ ), so every term with $j \geq 2$ is zero and the expansion has only two terms:

\mathbf{A}_{\text{II}}^k = 2^k \mathbf{I} + k \cdot 2^{k-1} \mathbf{N} = \begin{pmatrix} 2^k & k \cdot 2^{k-1} \\ 0 & 2^k \end{pmatrix}

(4)

The $k \cdot 2^{k-1}$ in the upper right corner is the source of “polynomial × exponential”—it is no accident; it is forced by the nilpotent part $\mathbf{N}$ .

	Form of $\mathbf{A}^k$	Structure of the solution
Matrix I ( $g=a=2$ )	$2^k \mathbf{I}$	Pure exponential: $\mathbf{x}_k = 2^k \mathbf{x}_0$
Matrix II ( $g<a$ )	$\begin{pmatrix} 2^k & k\cdot2^{k-1} \\ 0 & 2^k \end{pmatrix}$	Polynomial × exponential: a $k \cdot 2^{k-1}$ term appears

# ============================================================
# Part 1: numerical computation and visualization of the difference equation
# ============================================================
print_header('Difference equation: comparing the structure of A^k')

# --- Step 1: check the formula for A^k ---
print_step('1', 'Check A^k (k = 1, 2, 3, 5, 10)')
for k in [1, 2, 3, 5, 10]:
    Ak_diag   = np.linalg.matrix_power(A_diag, k)
    Ak_jordan = np.linalg.matrix_power(A_jordan, k)
    # analytic value from the Jordan block formula
    formula   = np.array([[2**k, k * 2**(k-1)],
                           [0,   2**k         ]], dtype=float)
    err = np.max(np.abs(Ak_jordan - formula))
    print(f'  k={k:2d}:  Matrix I A^k[0,1]={Ak_diag[0,1]:8.1f},  '
          f'Matrix II A^k[0,1]={Ak_jordan[0,1]:8.1f},  '
          f'formula error={err:.2e}')

# --- Step 2: trace the trajectory for the initial condition x0 = [1, 1] ---
print_step('2', 'Trace the trajectory: initial state x0 = [1, 1]')
x0 = np.array([1., 1.])
K_MAX = 15
k_vals = np.arange(K_MAX + 1)

# compute each step by iteration
traj_diag   = np.zeros((K_MAX + 1, 2))
traj_jordan = np.zeros((K_MAX + 1, 2))
traj_diag[0]   = x0
traj_jordan[0] = x0
for k in range(1, K_MAX + 1):
    traj_diag[k]   = A_diag   @ traj_diag[k-1]    # iterate step by step
    traj_jordan[k] = A_jordan @ traj_jordan[k-1]   # iterate step by step

# --- Step 3: compute the "growth rate": |x_k| / 2^k ---
print_step('3', 'Normalized growth: ||x_k|| / 2^k (pure exponential growth should approach a constant)')
norms_diag   = np.linalg.norm(traj_diag,   axis=1)
norms_jordan = np.linalg.norm(traj_jordan, axis=1)
growth_diag   = norms_diag   / (2.0**k_vals)
growth_jordan = norms_jordan / (2.0**k_vals)
print('  k:  ', '  '.join(f'{k:6d}' for k in k_vals[:8]))
print('  I:  ', '  '.join(f'{v:6.3f}' for v in growth_diag[:8]))
print('  II: ', '  '.join(f'{v:6.3f}' for v in growth_jordan[:8]))
print('  (I approaches a constant; the ratio for II grows linearly in k → faster than pure exponential growth)')

# --- Step 4: draw the comparison plots ---
print_step('4', 'Visualization: comparing the solutions of the difference equation')
fig = plt.figure(figsize=(14, 5), facecolor='white')
gs  = GridSpec(1, 3, figure=fig, wspace=0.38)

# Subplot 1: how the x_0 component changes with k
ax1 = fig.add_subplot(gs[0])
ax1.set_facecolor(LIGHT_BG)
ax1.plot(k_vals, traj_diag[:, 0],   'o-', color=C_VIOLET,
         lw=2, ms=5, label=r'I: $x_0$ component')
ax1.plot(k_vals, traj_jordan[:, 0], 's--', color=C_BLUE,
         lw=2, ms=5, label=r'II: $x_0$ component')
ax1.set_xlabel('Iteration step $k$', fontsize=11)
ax1.set_ylabel('Value of the $x_0$ component', fontsize=11)
ax1.set_title(r'$x_0$ component: I vs II', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Subplot 2: normalized growth ||x_k|| / 2^k
ax2 = fig.add_subplot(gs[1])
ax2.set_facecolor(LIGHT_BG)
ax2.plot(k_vals, growth_diag,   'o-', color=C_VIOLET,
         lw=2, ms=5, label=r'I: $\|\mathbf{x}_k\|/2^k$')
ax2.plot(k_vals, growth_jordan, 's--', color=C_BLUE,
         lw=2, ms=5, label=r'II: $\|\mathbf{x}_k\|/2^k$')
ax2.set_xlabel('Iteration step $k$', fontsize=11)
ax2.set_ylabel(r'$\|\mathbf{x}_k\|\,/\,2^k$', fontsize=11)
ax2.set_title(r'Normalized growth rate (divided by $2^k$)', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Subplot 3: growth of A^k[0,1] (the direct fingerprint of the Jordan structure)
ax3 = fig.add_subplot(gs[2])
ax3.set_facecolor(LIGHT_BG)
Ak_01_diag   = [np.linalg.matrix_power(A_diag,   k)[0,1] for k in k_vals]
Ak_01_jordan = [np.linalg.matrix_power(A_jordan, k)[0,1] for k in k_vals]
formula_vals = [k * 2**(k-1) if k > 0 else 0 for k in k_vals]   # formula k·2^{k-1}
ax3.plot(k_vals, Ak_01_diag,   'o-',  color=C_VIOLET, lw=2, ms=5,
         label=r'I: $\mathbf{A}^k_{01}=0$')
ax3.plot(k_vals, Ak_01_jordan, 's--', color=C_BLUE,   lw=2, ms=5,
         label=r'II: $\mathbf{A}^k_{01}=k\cdot2^{k-1}$')
ax3.plot(k_vals, formula_vals, 'x',   color=C_ORANGE, ms=8, lw=0,
         label='Formula check ✓')
ax3.set_xlabel('Iteration step $k$', fontsize=11)
ax3.set_ylabel(r'Upper right entry of $\mathbf{A}^k$', fontsize=11)
ax3.set_title(r'Fingerprint of the Jordan structure: $\mathbf{A}^k_{[0,1]}$', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

fig.suptitle('Solutions of the difference equation $\\mathbf{x}_{k+1} = \\mathbf{A}\\mathbf{x}_k$: diagonalizable vs Jordan block',
             fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

Part 2: The Differential Equation $\frac{d \mathbf{x} }{dt} = \mathbf{A}\mathbf{x}$ ¶

2.1 The Case of Matrix I: The Similarity Transformation Decouples Completely¶

For $\mathbf{A}_{\text{I}} = 2\mathbf{I}$ , the system is already decoupled:

\frac{d {x}_0 }{dt} = 2x_0, \qquad \frac{d {x}_1 }{dt} = 2x_1

(5)

The two directions are completely independent, and the solution is a pure exponential:

\mathbf{x}(t) = e^{2t}\mathbf{x}_0 = e^{2t}\begin{pmatrix}x_0(0) \\ x_1(0)\end{pmatrix}

(6)

This is precisely the geometric meaning of $g_2 = 2$ : there are two genuinely independent eigendirections, and the system evolves along each of them separately as $e^{2t}$ , without interference.

2.2 The Case of Matrix II: Incomplete Decoupling and the Appearance of $te^{2t}$ ¶

For $\mathbf{A}_{\text{II}} = \begin{pmatrix}2&1\\0&2\end{pmatrix}$ , the system is

\frac{d {x}_0 }{dt} = 2x_0 + x_1, \qquad \frac{d {x}_1 }{dt} = 2x_1

(7)

The equation for $x_1$ (the second one) can be solved on its own: $x_1(t) = x_1(0)\,e^{2t}$ . The first equation, however, contains $x_1$ , which must be substituted in before solving:

\frac{d {x}_0 }{dt} = 2x_0 + x_1(0)\,e^{2t}

(8)

By the method of variation of parameters, the solution is

x_0(t) = \bigl(x_0(0) + x_1(0)\,t\bigr)\,e^{2t}

(9)

(Note: differentiating the latter factor gives $2x_0$ , and differentiating the former factor gives $x_1(0)\,e^{2t}$ .)

Comparison of the complete solutions:

	$x_0(t)$	$x_1(t)$
Matrix I	$x_0(0)\,e^{2t}$	$x_1(0)\,e^{2t}$
Matrix II	$\bigl(x_0(0) + x_1(0)\,t\bigr)\,e^{2t}$	$x_1(0)\,e^{2t}$

# ============================================================
# Part 2: numerical computation and visualization of the differential equation
# ============================================================
print_header('Differential equation: analytic solution vs numerical solution by the matrix exponential')

# --- Initial condition and time axis ---
x0_ode = np.array([1., 1.])              # initial condition x(0) = [1, 1]
t_vals = np.linspace(0, 2.5, 300)       # time axis: 0 to 2.5

# --- Step 1: compute the matrix exponential solution with scipy.linalg.expm ---
print_step('1', 'Compute x(t) = e^{At} x0 (matrix exponential method)')
sol_diag   = np.array([scipy.linalg.expm(A_diag   * t) @ x0_ode for t in t_vals])
sol_jordan = np.array([scipy.linalg.expm(A_jordan * t) @ x0_ode for t in t_vals])

# --- Step 2: the analytic formulas ---
print_step('2', 'Compare with the analytic formulas')
x00, x10 = x0_ode
# analytic solution for Matrix I
x0_diag_exact = x00 * np.exp(2 * t_vals)
x1_diag_exact = x10 * np.exp(2 * t_vals)
# analytic solution for Matrix II
x0_jordan_exact = (x00 + x10 * t_vals) * np.exp(2 * t_vals)   # contains t·e^{2t}
x1_jordan_exact =  x10 * np.exp(2 * t_vals)

# check the error
err_diag   = np.max(np.abs(sol_diag[:,0]   - x0_diag_exact))
err_jordan = np.max(np.abs(sol_jordan[:,0] - x0_jordan_exact))
compare_print('Maximum error, analytic formula vs matrix exponential (x0 component):',
              f'{err_diag:.2e}', f'{err_jordan:.2e}')

# --- Step 3: normalized growth x(t) / e^{2t} ---
print_step('3', 'Normalized growth: x0(t) / e^{2t} (a pure exponential should give a constant)')
norm_diag_0   = sol_diag[:,0]   / np.exp(2 * t_vals)
norm_jordan_0 = sol_jordan[:,0] / np.exp(2 * t_vals)
print(f'  Matrix I: range of x0(t)/e^{{2t}}: [{norm_diag_0.min():.4f}, {norm_diag_0.max():.4f}]')
print(f'  Matrix II: range of x0(t)/e^{{2t}}: [{norm_jordan_0.min():.4f}, {norm_jordan_0.max():.4f}]')
print('  (I is the constant 1.0; II grows linearly from 1 as x0(0) + x1(0)·t)')

# --- Step 4: visualization ---
print_step('4', 'Visualization: comparing the solutions of the differential equation')
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5), facecolor='white')

# Subplot 1: direct comparison of x_0(t)
ax = axes[0]
ax.set_facecolor(LIGHT_BG)
ax.plot(t_vals, sol_diag[:,0],   color=C_VIOLET, lw=2.5, label='I: $x_0(0)e^{2t}$')
ax.plot(t_vals, sol_jordan[:,0], color=C_BLUE,   lw=2.5, ls='--',
        label='II: $(x_0(0)+x_1(0)t)e^{2t}$')
ax.set_xlabel('Time $t$', fontsize=11)
ax.set_ylabel('$x_0(t)$', fontsize=11)
ax.set_title('$x_0(t)$: I vs II', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Subplot 2: normalized growth x_0(t) / e^{2t}
ax = axes[1]
ax.set_facecolor(LIGHT_BG)
ax.plot(t_vals, norm_diag_0,   color=C_VIOLET, lw=2.5, label=r'I: constant')
ax.plot(t_vals, norm_jordan_0, color=C_BLUE,   lw=2.5, ls='--',
        label=r'II: $x_0(0)+x_1(0)\,t$ (linear growth)')
ax.set_xlabel('Time $t$', fontsize=11)
ax.set_ylabel(r'$x_0(t)\,/\,e^{2t}$', fontsize=11)
ax.set_title(r'Normalized growth (divided by $e^{2t}$)', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Subplot 3: phase-space trajectories (x0, x1)
ax = axes[2]
ax.set_facecolor(LIGHT_BG)
ax.plot(sol_diag[:,0],   sol_diag[:,1],   color=C_VIOLET, lw=2.5,
        label='I: straight-line trajectory')
ax.plot(sol_jordan[:,0], sol_jordan[:,1], color=C_BLUE,   lw=2.5, ls='--',
        label='II: curved trajectory')
ax.plot(*x0_ode, 'ko', ms=8, zorder=5, label='Initial point $\\mathbf{x}_0$')
ax.set_xlabel('$x_0$', fontsize=11)
ax.set_ylabel('$x_1$', fontsize=11)
ax.set_title('Phase-space trajectories $(x_0, x_1)$', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

fig.suptitle(r'Solutions of the differential equation $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}$: diagonalizable vs Jordan block',
             fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

Part 3: A Unifying View—The Matrix Exponential $e^{\mathbf{A}t}$ ¶

3.1 The Definition of the Matrix Exponential and Jordan Blocks¶

For any square matrix $\mathbf{A}$ , the matrix exponential is defined by the Taylor series

e^{\mathbf{A}t} = \sum_{k=0}^{\infty} \frac{(\mathbf{A}t)^k}{k!} = \mathbf{I} + \mathbf{A}t + \frac{\mathbf{A}^2 t^2}{2!} + \cdots

(10)

Diagonal matrix: $e^{2\mathbf{I}\cdot t} = e^{2t}\mathbf{I}$ , a pure exponential; every term is a multiple of the identity matrix.

Jordan block $\mathbf{A}_{\text{II}} = 2\mathbf{I} + \mathbf{N}$ , $\mathbf{N}^2 = \mathbf{0}$ :

e^{(2\mathbf{I}+\mathbf{N})t} = e^{2\mathbf{I}t} \cdot e^{\mathbf{N}t} = e^{2t}\mathbf{I} \cdot \left(\mathbf{I} + \mathbf{N}t + \underbrace{\frac{\mathbf{N}^2 t^2}{2!}}_{=\,\mathbf{0}} + \cdots\right) = e^{2t}\begin{pmatrix}1 & t \\ 0 & 1\end{pmatrix}

(11)

3.2 Unifying the Difference Equation and the Differential Equation¶

Case	Matrix exponential $e^{\mathbf{A}t}$	Difference equation $\mathbf{A}^k$	Reason
$\mathbf{A}_{\text{I}} = 2\mathbf{I}$	$e^{2t}\mathbf{I}$	$2^k \mathbf{I}$	$\mathbf{N}=\mathbf{0}$ , no extra structure
$\mathbf{A}_{\text{II}} = 2\mathbf{I}+\mathbf{N}$	$e^{2t}\begin{pmatrix}1&t\\0&1\end{pmatrix}$	$\begin{pmatrix}2^k & k\cdot2^{k-1}\\0&2^k\end{pmatrix}$	$\mathbf{N}^2=\mathbf{0}$ , the expansion terminates at the first-degree term

# ============================================================
# Part 3: checking and visualizing the matrix exponential
# ============================================================
print_header('Comparing the structure of the matrix exponential e^{At}')

# --- Step 1: check the matrix exponential formula ---
print_step('1', 'Check the analytic formula for e^{At} against scipy.linalg.expm')
t_check = [0.0, 0.5, 1.0, 2.0]
print('  {:>6s}  {:>20s}  {:>20s}  {:>12s}'.format(
      't', 'I expm[0,1]', 'II expm[0,1]', 'II formula t·e^{2t}'))
for t in t_check:
    expm_d = scipy.linalg.expm(A_diag   * t)
    expm_j = scipy.linalg.expm(A_jordan * t)
    formula_j01 = t * np.exp(2 * t)           # formula: e^{2t} times the t in the upper right corner
    print(f'  {t:6.2f}  {expm_d[0,1]:20.8f}  {expm_j[0,1]:20.8f}  {formula_j01:12.8f}')

# --- Step 2: visualizing the termination of the Taylor expansion ---
print_step('2', 'Termination of the Taylor expansion: the contribution of each term of e^{Nt}')
t0 = 1.5    # fixed time
N  = np.array([[0., 1.], [0., 0.]])           # nilpotent part
print(f'  N   = {N[0]}')
print(f'        {N[1]}')
print(f'  N^2 = {(N@N)[0]}')
print(f'        {(N@N)[1]}')
print(f'  t={t0}: I + N·t = I + {t0}·N')
term0 = np.eye(2)
term1 = N * t0
term2 = (N @ N) * t0**2 / 2          # = 0
print(f'  Term 0 (I):\n    {term0}')
print(f'  Term 1 (N·t):\n    {term1}')
print(f'  Term 2 (N²t²/2):\n    {term2}  ← always zero!')
e_Nt = term0 + term1                  # only the first two terms are needed
e_At_jordan = np.exp(2*t0) * e_Nt   # the full e^{At}
compare_print(
    f'e^{{At}} at t={t0} (first two terms only vs scipy):',
    f'\n    {scipy.linalg.expm(A_diag*t0).round(6)}',
    f'\n    formula={e_At_jordan.round(6)}'  
    f'\n    scipy={scipy.linalg.expm(A_jordan*t0).round(6)}')

# --- Step 3: plot how each entry of the matrix exponential evolves in time ---
print_step('3', 'Visualization: time evolution of the entries of the matrix exponential e^{At}')
t_vals2 = np.linspace(0, 2.5, 200)

# compute the entries of the matrix exponential at all times
expm_diag_00   = np.exp(2 * t_vals2)          # [0,0] entry
expm_diag_01   = np.zeros_like(t_vals2)        # [0,1] entry (always 0)
expm_jordan_00 = np.exp(2 * t_vals2)           # [0,0] entry (the same)
expm_jordan_01 = t_vals2 * np.exp(2 * t_vals2) # [0,1] entry (= t·e^{2t})

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), facecolor='white')

# Subplot 1: the [0,0] entry (the same for both)
ax = axes[0]
ax.set_facecolor(LIGHT_BG)
ax.plot(t_vals2, expm_diag_00, color=C_VIOLET, lw=3,
        label=r'I / II: $e^{\mathbf{A}t}_{[0,0]} = e^{2t}$')
ax.set_xlabel('Time $t$', fontsize=11)
ax.set_ylabel(r'$e^{\mathbf{A}t}_{[0,0]}$', fontsize=11)
ax.set_title(r'Diagonal entry (the same for both matrices)', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Subplot 2: the [0,1] entry (the key difference)
ax = axes[1]
ax.set_facecolor(LIGHT_BG)
ax.axhline(0, color=C_VIOLET, lw=3,
           label=r'I: $e^{\mathbf{A}t}_{[0,1]} = 0$ (no Jordan structure)')
ax.plot(t_vals2, expm_jordan_01, color=C_BLUE, lw=3, ls='--',
        label=r'II: $e^{\mathbf{A}t}_{[0,1]} = t\,e^{2t}$ (Jordan structure)')

# mark the contribution of the Taylor expansion
ax.fill_between(t_vals2, 0, expm_jordan_01,
                alpha=0.15, color=C_BLUE, label='contribution of $t\\cdot e^{2t}$')
ax.set_xlabel('Time $t$', fontsize=11)
ax.set_ylabel(r'$e^{\mathbf{A}t}_{[0,1]}$', fontsize=11)
ax.set_title(r'Upper right entry: the fingerprint of the Jordan structure', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

fig.suptitle(r'Comparing the entries of the matrix exponential $e^{\mathbf{A}t}$: $\mathbf{N}^2=\mathbf{0}$ makes the expansion terminate',
             fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

3.3 A Comparison of the Three Views¶

Unified Summary

For the same “repeated eigenvalue $\lambda=2$ with algebraic multiplicity $=2$ ,” the gap in geometric multiplicity shows up in exactly the same way in three different views:

View	Diagonalizable ( $g=a=2$ )	Jordan block ( $g<a$ )
Eigenspace	2 independent eigendirections	Only 1; one generalized eigenvector is needed to make up the difference
Structure of $\mathbf{N}$	$\mathbf{N}=\mathbf{0}$	$\mathbf{N}\neq\mathbf{0}$ , $\mathbf{N}^2=\mathbf{0}$
Difference equation $\mathbf{A}^k$	$2^k\mathbf{I}$ (pure exponential)	$k\cdot2^{k-1}$ appears (polynomial × exponential)
Differential equation $e^{\mathbf{A}t}$	$e^{2t}\mathbf{I}$ (pure exponential)	$t\cdot e^{2t}$ appears (polynomial × exponential)
Long-term behavior	Grows at the rate $2^k$ (or $e^{2t}$ )	Grows by an extra factor of $k$ (or $t$ )

“Polynomial × exponential” is not an accidental outcome of the computation; it is the inevitable manifestation of the termination after finitely many terms of the Taylor expansion when $\mathbf{N}^2=\mathbf{0}$ .

Observation and Reflection: The Matrix Exponential

In Step 2, the Taylor expansion of $e^{\mathbf{N}t}$ terminates at the second term because $\mathbf{N}^2=\mathbf{0}$ . If the Jordan block were of size $3\times3$ (that is, $\mathbf{N}^3=\mathbf{0}$ but $\mathbf{N}^2\neq\mathbf{0}$ ), at which term would the expansion terminate? Which types of terms would appear in the solution?
For Matrix II, why can $e^{(2\mathbf{I}+\mathbf{N})t}$ be written as $e^{2\mathbf{I}t}\cdot e^{\mathbf{N}t}$ ? What condition does this require? (Hint: check that $2\mathbf{I}\cdot\mathbf{N} = \mathbf{N}\cdot 2\mathbf{I}$ .)
According to the table in the unified summary, the “ $k$ ” in the difference equation and the “ $t$ ” in the differential equation play the same role. What is their underlying source?

Part 4: Exercises¶

The following exercises are set in the context of radioactive decay. They form three progressive problems, from the concrete to the general, echoing the decay chains in §8.3.3 and §8.4.2.

Exercise 1 (Different Decay Rates: Decoupling by a Similarity Transformation)

Two isotopes A and B form the decay chain A → B → C (C is stable), with decay rate constants $k_A = 3$ and $k_B = 1$ :

\frac{d}{dt}\begin{pmatrix}x_0 \\ x_1\end{pmatrix} = \underbrace{\begin{pmatrix}-3 & 0 \\ 3 & -1\end{pmatrix}}_{\mathbf{A}} \begin{pmatrix}x_0 \\ x_1\end{pmatrix}, \qquad \mathbf{x}(0) = \begin{pmatrix}1 \\ 0\end{pmatrix}

(12)

Find the eigenvalues and eigenvectors of $\mathbf{A}$ , and confirm that $\mathbf{A}$ is diagonalizable.
Use the similarity transformation $\mathbf{y} = \mathbf{P}^{-1}\mathbf{x}$ to decouple the system, and write down the solution of $\dot{\mathbf{y}} = \mathbf{D}\mathbf{y}$ .
Substitute back to obtain explicit expressions for $x_0(t)$ and $x_1(t)$ .
Verify the result in Python with scipy.linalg.expm, and plot $x_0(t)$ and $x_1(t)$ as functions of time.
Find the time $t^*$ at which $x_1(t)$ attains its maximum.

Exercise 2 (Equal Decay Rates: A Jordan Block and

te^{-t}

)

Now let $k_A = k_B = 1$ (the two decay rates are equal):

\frac{d}{dt}\begin{pmatrix}x_0 \\ x_1\end{pmatrix} = \underbrace{\begin{pmatrix}-1 & 0 \\ 1 & -1\end{pmatrix}}_{\mathbf{B}} \begin{pmatrix}x_0 \\ x_1\end{pmatrix}, \qquad \mathbf{x}(0) = \begin{pmatrix}1 \\ 0\end{pmatrix}

(16)

Find the eigenvalue of $\mathbf{B}$ together with its algebraic and geometric multiplicities, and confirm that $\mathbf{B}$ is not diagonalizable.
Use $\mathbf{B} = -\mathbf{I} + \mathbf{M}$ (where $\mathbf{M} = \begin{pmatrix}0&0\\1&0\end{pmatrix}$ ) and the binomial expansion to compute $e^{\mathbf{B}t}$ directly.
Solve for $x_0(t)$ and $x_1(t)$ , and verify that $x_1(t)$ contains a $te^{-t}$ term.
Compare with the result of Exercise 1: with equal decay rates, does the time $t^{**}$ at which $x_1(t)$ attains its maximum change? (Hint: differentiate $te^{-t}$ .)
In Python, plot the $x_1(t)$ of Exercise 1 and of Exercise 2 on the same figure to compare the two cases visually.

Exercise 3 (◆ Exploration: Geometric Multiplicity and the Variety of Solutions)

Consider a family of matrices (with parameter $\epsilon \geq 0$ ):

\mathbf{C}(\epsilon) = \begin{pmatrix}2 & \epsilon \\ 0 & 2\end{pmatrix}

(18)

For every $\epsilon$ , the algebraic multiplicity of $\mathbf{C}(\epsilon)$ is 2. When is the geometric multiplicity 1, and when is it 2?
For the initial condition $\mathbf{x}(0)=[1,1]^\top$ and the difference equation $\mathbf{x}_{k+1}=\mathbf{C}(\epsilon)\mathbf{x}_k$ , what is the closed-form expression for the $x_0$ component of $\mathbf{x}_k$ ? (Express it in terms of $\epsilon$ and $k$ .)
Take $\epsilon = 0, 0.1, 0.5, 1.0, 2.0$ and use Python to plot the trajectories of the $x_0$ component for $k=0,1,\ldots,12$ (on the same figure). Observe how, as $\epsilon$ increases, the trajectory passes from “pure exponential” to “polynomial × exponential.”
Explain: why does the effect of the Jordan structure disappear as $\epsilon \to 0$ , whereas for $\epsilon > 0$ , no matter how small, the long-term growth outpaces the pure exponential $2^k$ ?

Key Takeaways of This Experiment

Geometric multiplicity = the number of independent decay modes: $g_\lambda$ independent eigendirections → $g_\lambda$ pure exponential modes; the “missing directions” are replaced by Jordan chains, which bring in polynomial factors.
Similarity transformation = changing language so that the equations decouple: $\dot{\mathbf{y}} = \mathbf{D}\mathbf{y}$ looks like a miracle of decoupling, but it is no miracle—it is $g_\lambda = a_\lambda$ that guarantees enough independent eigendirections.
The source of $te^{\lambda t}$ : the Taylor expansion $e^{\mathbf{N}t} = \mathbf{I} + \mathbf{N}t$ of the nilpotent matrix $\mathbf{N}$ (which terminates when $\mathbf{N}^2=\mathbf{0}$ ), multiplied by the scalar exponential $e^{\lambda t}$ , naturally produces $te^{\lambda t}$ .
The three views agree: the difference between the algebraic and geometric multiplicities gives exactly the same information about “the degree of the polynomial factor” in all three views—eigenspaces, difference equations, and differential equations.