14 — Control Systems Stability: Complex Eigenvalues & EIGEN_GENERAL

Two real linear-time-invariant (LTI) systems, both genuinely non-symmetric — the routine case for real aircraft, circuits, and mechanical systems, which EIGENVALUES/EIGEN (symmetric-only) simply cannot handle at all. This notebook is the deferred validation notebook for Phase 3 of linal-db-rs/SCIENTIFIC_ENGINE_EXPANSION_PLAN.md: the scalar Value::Complex type (COMPLEX(re, im), REAL/IMAG/ABS/PHASE/CONJ) and EIGENVALUES_GENERAL/EIGEN_GENERAL (general, possibly-non-symmetric eigendecomposition via nalgebra's Schur decomposition), released in engine v0.1.84 / linaldb 0.1.11.

For an LTI system dx/dt = A x, an eigenvalue's real part is the decay/growth rate (negative = stable) and its imaginary part is the oscillation angular frequency (zero = non-oscillatory, purely exponential decay/growth). A real, non-symmetric matrix can have either real or complex (conjugate-pair) eigenvalues — you don't know which until you compute them, and the previous symmetric-only EIGEN/EIGENVALUES couldn't even attempt a general matrix.

Like notebook 11 (11_dynamical_systems_chaos.ipynb), this notebook uses real, citable physics/mathematics rather than a downloaded dataset — a real published aircraft flight-dynamics model and a real circuit built from standard component values, both cross-checked against independent numpy computations and, for the aircraft, against the original published source's own numbers.

import cmath
import math

import numpy as np

import linaldb

print("linaldb version:", linaldb.__version__)

db = linaldb.Db(data_dir="./data_control_systems_stability")

linaldb version: 0.1.12

Part 1 — Boeing 747 longitudinal flight dynamics (complex eigenvalues)

Source: Stephen Boyd, EE263 (Stanford), Lecture 14: "Example: Aircraft dynamics" — a real, standard, widely-cited control-theory teaching example with an exact published numeric state matrix. For a Boeing 747 in level flight at 40,000 ft, 774 ft/sec:

$$\dot{x} = Ax, \quad x = [u, v, q, \theta]$$

where u = forward-velocity deviation (ft/s), v = perpendicular-velocity deviation (ft/s, down positive), q = pitch rate (crad/s), theta = pitch angle (crad, up positive) — all small deviations from trim. crad = centiradian = 0.01 rad ≈ 0.57°.

$$A = \begin{bmatrix} -0.003 & 0.039 & 0 & -0.322 \\ -0.065 & -0.319 & 7.74 & 0 \\ 0.020 & -0.101 & -0.429 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}$$

This matrix is not symmetric (e.g. A[0][3] = -0.322 vs A[3][0] = 0), so plain EIGENVALUES/EIGEN cannot touch it at all.

db.execute(
    "MATRIX boeing747 = [[-0.003, 0.039, 0, -0.322], "
    "[-0.065, -0.319, 7.74, 0], [0.020, -0.101, -0.429, 0], [0, 0, 1, 0]]"
)
db.execute("LET eigs747 = EIGENVALUES_GENERAL boeing747")
eigs747 = db.execute("SHOW eigs747").to_numpy()
print("linaldb EIGENVALUES_GENERAL -- shape:", eigs747.shape)
print("row 0 (real parts):     ", eigs747[0])
print("row 1 (imaginary parts):", eigs747[1])

linaldb EIGENVALUES_GENERAL -- shape: (2, 4)
row 0 (real parts):      [-0.37504214 -0.37504214 -0.00045786 -0.00045786]
row 1 (imaginary parts): [ 0.88175201 -0.88175201  0.06737732 -0.06737732]

linaldb_eigs = sorted(
    (complex(re, im) for re, im in zip(eigs747[0], eigs747[1])),
    key=lambda z: (z.real, z.imag),
)

A747 = np.array(
    [
        [-0.003, 0.039, 0, -0.322],
        [-0.065, -0.319, 7.74, 0],
        [0.020, -0.101, -0.429, 0],
        [0, 0, 1, 0],
    ]
)
numpy_eigs = sorted(np.linalg.eigvals(A747), key=lambda z: (z.real, z.imag))

# Boyd's own published values (EE263 Lecture 14, slide "Eigenvalues and modes"):
# -0.3750 +/- 0.8818j (short-period), -0.0005 +/- 0.0674j (phugoid)
boyd_published = sorted(
    [
        complex(-0.3750, -0.8818), complex(-0.3750, 0.8818),
        complex(-0.0005, -0.0674), complex(-0.0005, 0.0674),
    ],
    key=lambda z: (z.real, z.imag),
)

print(f"{'linaldb EIGENVALUES_GENERAL':32s} {'numpy eigvals':32s} {'Boyd (EE263, published)':24s}")
for a, b, c in zip(linaldb_eigs, numpy_eigs, boyd_published):
    print(f"{a:+.6f}{'':6s}      {b:+.6f}{'':6s}      {c:+.4f}")
    assert abs(a - b) < 1e-4
    assert abs(a - c) < 5e-3

print("\nAll three sources agree.")

linaldb EIGENVALUES_GENERAL      numpy eigvals                    Boyd (EE263, published) 
-0.375042-0.881752j            -0.375042-0.881752j            -0.3750-0.8818j
-0.375042+0.881752j            -0.375042+0.881752j            -0.3750+0.8818j
-0.000458-0.067377j            -0.000458-0.067377j            -0.0005-0.0674j
-0.000458+0.067377j            -0.000458+0.067377j            -0.0005+0.0674j

All three sources agree.

Physical interpretation: short-period and phugoid modes

Both complex-conjugate pairs have negative real parts (stable) but nonzero imaginary parts (oscillatory) — genuine, well-known aircraft handling-qualities concepts. Real part = decay rate, imaginary part = angular frequency, so period = 2*pi / |imag| and the "1/e" decay time = 1 / |real|.

We compute these directly from the complex eigenvalues using the engine's own Complex scalar type — COMPLEX(re, im) plus ABS/PHASE — not just consuming EIGENVALUES_GENERAL's raw Matrix(2,N) output.

db.execute(
    "DATASET eigs747_tbl COLUMNS (mode: String, re: Double, im: Double)"
)
short_re, short_im = float(eigs747[0][0]), float(abs(eigs747[1][0]))
phug_idx = int(np.argmin(np.abs(eigs747[0])))
phug_re, phug_im = float(eigs747[0][phug_idx]), float(abs(eigs747[1][phug_idx]))

db.execute(
    f"INSERT INTO eigs747_tbl (mode = 'short_period', re = {short_re!r}, im = {short_im!r})"
)
db.execute(
    f"INSERT INTO eigs747_tbl (mode = 'phugoid', re = {phug_re!r}, im = {phug_im!r})"
)

result = db.execute(
    "SELECT mode, re AS decay_rate, "
    "ABS(COMPLEX(re, im)) AS natural_freq_rad_s, "
    "PHASE(COMPLEX(re, im)) AS phase_rad "
    "FROM eigs747_tbl"
)
print(f"{'mode':14s} {'decay rate (1/s)':18s} {'|z| (rad/s)':14s} {'period (s)':12s} {'1/e decay (s)':14s}")
for mode, decay_rate, wn, phase in result.rows:
    imag = wn * math.sin(phase)  # reconstruct from |z| and phase = atan2(im, re)
    period = 2 * math.pi / abs(imag)
    decay_time = 1 / abs(decay_rate)
    print(f"{mode:14s} {decay_rate:+18.6f} {wn:14.6f} {period:12.2f} {decay_time:14.2f}")

    # cross-check ABS/PHASE against Python's own cmath on the identical complex number
    z = complex(decay_rate, imag)
    assert abs(wn - abs(z)) < 1e-6
    assert abs(phase - cmath.phase(z)) < 1e-6

print("\nlinaldb's ABS/PHASE match Python's cmath exactly on the same complex values.")
print("Short-period: ~7s period, decays in a few seconds -- barely noticeable to a pilot.")
print("Phugoid: ~90s period, decays over thousands of seconds -- a slow, gentle")
print("  speed/altitude oscillation pilots learn to recognize and damp out manually.")

mode           decay rate (1/s)   |z| (rad/s)    period (s)   1/e decay (s) 
short_period            -0.375042       0.958198         7.13           2.67
phugoid                 -0.000458       0.067379        93.25        2184.05

linaldb's ABS/PHASE match Python's cmath exactly on the same complex values.
Short-period: ~7s period, decays in a few seconds -- barely noticeable to a pilot.
Phugoid: ~90s period, decays over thousands of seconds -- a slow, gentle
  speed/altitude oscillation pilots learn to recognize and damp out manually.

EIGEN_GENERAL correctly refuses this matrix

EIGEN_GENERAL only computes real eigenvectors (nalgebra has no public general complex-eigenvector solver). Since this matrix's eigenvalues are genuinely complex, it must error loudly here rather than silently returning wrong or partial vectors — exactly the honest-constraint pattern notebook 11 documents for EIGEN's symmetric-only restriction.

try:
    db.execute("LET vals747, vecs747 = EIGEN_GENERAL boeing747")
    raise AssertionError("expected EIGEN_GENERAL to reject genuinely complex eigenvalues")
except linaldb.LinalError as e:
    print(f"EIGEN_GENERAL on the aircraft matrix correctly errors:\n  {e}")

EIGEN_GENERAL on the aircraft matrix correctly errors:
  [line 8] Engine error: Invalid operation: EIGEN_GENERAL: matrix has complex eigenvalues -- eigenvector computation for the complex case is not implemented; use EIGENVALUES_GENERAL for the eigenvalues alone

Part 2 — A two-stage RC lowpass filter cascade (real eigenvalues, non-symmetric)

Real physics (Kirchhoff's current law), real standard E12-series component values — deliberately different per stage (R1=10kΩ, C1=100nF; R2=4.7kΩ, C2=220nF), unlike a matched-stage cascade, which would give a symmetric matrix and wouldn't need EIGEN_GENERAL at all.

State [v1, v2] = voltage at the node between the two stages, voltage at the output node (after stage 2), both relative to ground, with input vin driving stage 1 through R1. Kirchhoff's current law at each node:

$$C_1 \dot{v}_1 = \frac{v_{in}-v_1}{R_1} - \frac{v_1-v_2}{R_2} \;\Rightarrow\; \dot{v}_1 = -\left(\frac{1}{R_1C_1}+\frac{1}{R_2C_1}\right)v_1 + \frac{1}{R_2C_1}v_2 + \frac{1}{R_1C_1}v_{in}$$

$$C_2 \dot{v}_2 = \frac{v_1-v_2}{R_2} \;\Rightarrow\; \dot{v}_2 = \frac{1}{R_2C_2}v_1 - \frac{1}{R_2C_2}v_2$$

giving the (input-free, homogeneous) system matrix A below.

R1, C1 = 10e3, 100e-9   # stage 1: 10 kOhm, 100 nF
R2, C2 = 4.7e3, 220e-9  # stage 2: 4.7 kOhm, 220 nF

a00 = -(1 / (R1 * C1) + 1 / (R2 * C1))
a01 = 1 / (R2 * C1)
a10 = 1 / (R2 * C2)
a11 = -1 / (R2 * C2)
print(f"A = [[{a00:.10f}, {a01:.10f}], [{a10:.10f}, {a11:.10f}]]")
print(f"symmetric? A[0][1]={a01:.6f} vs A[1][0]={a10:.6f} -- no, genuinely different")

db.execute(f"MATRIX rcfilter = [[{a00!r}, {a01!r}], [{a10!r}, {a11!r}]]")

A = [[-3127.6595744681, 2127.6595744681], [967.1179883946, -967.1179883946]]
symmetric? A[0][1]=2127.659574 vs A[1][0]=967.117988 -- no, genuinely different

'Defined matrix: rcfilter'

db.execute("LET rcvals, rcvecs = EIGEN_GENERAL rcfilter")
rcvals = db.execute("SHOW rcvals").to_numpy()
rcvecs = db.execute("SHOW rcvecs").to_numpy()
print("linaldb EIGEN_GENERAL eigenvalues: ", rcvals)
print("linaldb EIGEN_GENERAL eigenvectors (columns):\n", rcvecs)

A_rc = np.array([[a00, a01], [a10, a11]])
np_vals, np_vecs = np.linalg.eig(A_rc)
print("\nnumpy eigenvalues: ", np_vals.real)
print("numpy eigenvectors (columns):\n", np_vecs.real)

# Match each linaldb eigenvalue to its closest numpy counterpart (order isn't
# guaranteed to agree between the two implementations) before comparing.
for i in range(2):
    j = int(np.argmin(np.abs(np_vals.real - rcvals[i])))
    assert abs(rcvals[i] - np_vals[j].real) < 1e-2
    # an eigenvector is only defined up to sign/scale -- compare direction
    cos_sim = abs(np.dot(rcvecs[:, i], np_vecs[:, j].real))
    assert cos_sim > 0.9999

print("\nBoth eigenvalues and eigenvector directions match numpy.")
print("Both eigenvalues are real and negative -- no resonance is physically possible")
print("in a passive RC network with no inductance, so this had to come out real.")

linaldb EIGEN_GENERAL eigenvalues:  [-3843.12915039  -251.64857483]
linaldb EIGEN_GENERAL eigenvectors (columns):
 [[-0.94784474 -0.5947367 ]
 [ 0.31873232 -0.80392057]]

numpy eigenvalues:  [-3843.12897253  -251.64859033]
numpy eigenvectors (columns):
 [[-0.94784476 -0.5947367 ]
 [ 0.31873234 -0.80392056]]

Both eigenvalues and eigenvector directions match numpy.
Both eigenvalues are real and negative -- no resonance is physically possible
in a passive RC network with no inductance, so this had to come out real.

Physical interpretation: two real time constants

$1/|\lambda|$ gives each mode's decay time constant — the fast and slow decay modes of the cascade. We also verify A v = \lambda v directly, both through the engine's own MATMUL and through numpy, as a real correctness check on the returned eigenvectors (not just their eigenvalues).

for i in range(2):
    tau = 1 / abs(rcvals[i])
    unit = "us" if tau < 1e-3 else "ms"
    tau_disp = tau * 1e6 if unit == "us" else tau * 1e3
    print(f"lambda_{i} = {rcvals[i]:+.4f} /s  ->  time constant = {tau_disp:.3f} {unit}")

# A v = lambda v, verified via linaldb's own MATMUL (matrix * column-vector-as-matrix)
v0 = rcvecs[:, 0]
db.execute(f"MATRIX v0col = [[{float(v0[0])!r}], [{float(v0[1])!r}]]")
db.execute("LET av0 = MATMUL rcfilter v0col")
av0 = db.execute("SHOW av0").to_numpy().ravel()
lv0 = rcvals[0] * v0

print(f"\nA @ v0 (linaldb MATMUL) = {av0}")
print(f"lambda0 * v0            = {lv0}")
assert np.allclose(av0, lv0, atol=1e-2)

# cross-check the same identity directly in numpy too
assert np.allclose(A_rc @ v0, lv0, atol=1e-6)
print("\nA @ v = lambda * v holds, both via linaldb's MATMUL and via numpy.")

lambda_0 = -3843.1292 /s  ->  time constant = 260.205 us
lambda_1 = -251.6486 /s  ->  time constant = 3.974 ms

A @ v0 (linaldb MATMUL) = [ 3642.68994141 -1224.92944336]
lambda0 * v0            = [ 3642.68976467 -1224.92947502]

A @ v = lambda * v holds, both via linaldb's MATMUL and via numpy.

Plain EIGEN correctly refuses this matrix

This is the direct motivating contrast for EIGEN_GENERAL: a real, physically legitimate system whose state matrix simply isn't symmetric.

try:
    db.execute("LET ev, vv = EIGEN rcfilter")
    raise AssertionError("expected EIGEN to reject a non-symmetric matrix")
except linaldb.LinalError as e:
    print(f"EIGEN on the RC filter matrix correctly errors:\n  {e}")

EIGEN on the RC filter matrix correctly errors:
  [line 16] Engine error: Invalid operation: EIGEN: matrix is not symmetric (entries [0][1]=2127.65966796875 vs [1][0]=967.1179809570313 differ) -- only symmetric matrices are supported today

Summary

• Value::Complex (COMPLEX(re, im), REAL/IMAG/ABS/PHASE/CONJ) exercised directly on real eigenvalue data — ABS/PHASE matched Python's cmath exactly on the same complex numbers.
• EIGENVALUES_GENERAL on the real Boeing 747 longitudinal-dynamics matrix (Stephen Boyd, EE263 Lecture 14) reproduced the published complex eigenvalues for both the short-period and phugoid modes, matching both the original published numbers and an independent numpy computation.
• EIGEN_GENERAL on that same matrix correctly refused to compute eigenvectors (genuinely complex eigenvalues, no general complex-eigenvector solver) — a loud, documented error, not a silently wrong answer.
• EIGEN_GENERAL on a real, non-symmetric two-stage RC filter cascade (derived from Kirchhoff's current law, standard E12 component values) correctly returned two real eigenvalues and eigenvectors, matching numpy in both value and direction, with A v = \lambda v verified directly via linaldb's own MATMUL.
• Plain EIGEN correctly refused that same non-symmetric RC-filter matrix — the direct motivating contrast for why EIGEN_GENERAL exists at all: real aircraft, circuits, and mechanical systems are routinely non-symmetric, and the previous symmetric-only restriction couldn't attempt them at all, whether their eigenvalues turned out real or complex.

Every numerical result in this notebook was cross-checked against either the original published source (Boyd's EE263 lecture) or a direct numpy computation — no cherry-picking.