Solving ODEs with Matrix Exponentials¶

The linear ODE $\dot{\mathbf{x}} = A\mathbf{x}$ has the exact solution $\mathbf{x}(t) = e^{At}\mathbf{x}_0$. This notebook computes matrix exponentials numerically with np_expm and compares them against symbolic solutions, building intuition for how eigenstructure governs solution behaviour.

In [181]:
load("numerics")$
load("ax-plots")$

The Method¶

Given a constant matrix $A$ and initial condition $\mathbf{x}_0$, the solution to

$$\frac{d\mathbf{x}}{dt} = A\mathbf{x}, \quad \mathbf{x}(0) = \mathbf{x}_0$$

is the matrix exponential:

$$\mathbf{x}(t) = e^{At}\,\mathbf{x}_0$$

The matrix exponential $e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots$ is computed numerically by np_expm (wrapping LAPACK's scaling-and-squaring algorithm). For each time step $t$, we form $At$, compute $e^{At}$, and multiply by $\mathbf{x}_0$.

Simple 2D Decoupled System [S+N]¶

We start with a diagonal (decoupled) system:

$$A = \begin{pmatrix} -1 & 0 \\ 0 & -2 \end{pmatrix}$$

The symbolic solution is $x(t) = x_0 e^{-t}$, $y(t) = y_0 e^{-2t}$. Both components decay exponentially, with $y$ decaying twice as fast.

In [182]:
/* Symbolic solution for the decoupled system */
A_sym : matrix([-1, 0], [0, -2])$
print("A =", A_sym)$

/* Eigenvalues confirm decoupled exponential decay */
evals : eigenvalues(A_sym);
print("Eigenvalues:", evals[1])$

/* Symbolic solution */
print("x(t) = x0 * exp(-t)")$
print("y(t) = y0 * exp(-2t)")$

A = matrix([-1,0],[0,-2])
[[-2,-1],[1,1]]
Eigenvalues: [-2,-1]
x(t) = x0 * exp(-t)
y(t) = y0 * exp(-2t)
In [183]:
/* Numeric trajectory via matrix exponential */
A_nd : ndarray(float(A_sym))$
x0 : ndarray([1.0, 1.0])$

ts : makelist(i*0.05, i, 0, 100)$
trajectories : makelist(
  block([M : np_expm(np_scale(t, A_nd))],
    np_to_list(np_matmul(M, x0))),
  t, ts)$

xs : makelist(first(p), p, trajectories)$
ys : makelist(second(p), p, trajectories)$

ax_draw2d(
  color="blue", lines(ts, xs), name="x(t)",
  color="red", lines(ts, ys), name="y(t)",
  title="Decoupled System: x(t) and y(t)",
  xlabel="t", ylabel="value",
  showlegend=true
)$
No description has been provided for this image
In [184]:
/* Phase portrait: curve in the (x, y) plane */
ax_draw2d(
  lines(xs, ys),
  title="Decoupled System: Phase Portrait",
  xlabel="x", ylabel="y"
)$
No description has been provided for this image

Spiral System¶

A coupled system with complex eigenvalues produces spiralling trajectories:

$$A = \begin{pmatrix} -0.5 & 1 \\ -1 & -0.5 \end{pmatrix}$$

The eigenvalues are $-0.5 \pm i$, giving damped oscillation: the solution spirals inward toward the origin.

In [185]:
/* Spiral system: damped rotation */
A_spiral : matrix([-0.5, 1], [-1, -0.5])$
A_spiral_nd : ndarray(float(A_spiral))$
x0_spiral : ndarray([2.0, 0.0])$

ts_spiral : makelist(i*0.05, i, 0, 200)$
traj_spiral : makelist(
  block([M : np_expm(np_scale(t, A_spiral_nd))],
    np_to_list(np_matmul(M, x0_spiral))),
  t, ts_spiral)$

xs_sp : makelist(first(p), p, traj_spiral)$
ys_sp : makelist(second(p), p, traj_spiral)$

ax_draw2d(
  lines(xs_sp, ys_sp),
  title="Damped Spiral: Phase Portrait",
  xlabel="x", ylabel="y",
  aspect_ratio=true
)$
No description has been provided for this image
In [186]:
/* Time series: oscillating decay */
ax_draw2d(
  color="blue", lines(ts_spiral, xs_sp), name="x(t)",
  color="red", lines(ts_spiral, ys_sp), name="y(t)",
  title="Damped Spiral: Time Series",
  xlabel="t", ylabel="value",
  showlegend=true
)$
No description has been provided for this image

3D Trajectory¶

A 3D system with rotation in the $(x,y)$ plane and independent decay in $z$:

$$A = \begin{pmatrix} -0.1 & 1 & 0 \\ -1 & -0.1 & 0 \\ 0 & 0 & -0.5 \end{pmatrix}$$

The trajectory spirals inward in $(x,y)$ while $z$ decays independently.

In [187]:
/* 3D spiral-decay system */
A_3d : matrix([-0.1, 1, 0], [-1, -0.1, 0], [0, 0, -0.5])$
A_3d_nd : ndarray(float(A_3d))$
x0_3d : ndarray([1.0, 0.0, 2.0])$

ts_3d : makelist(i*0.1, i, 0, 200)$
traj_3d : makelist(
  block([M : np_expm(np_scale(t, A_3d_nd))],
    np_to_list(np_matmul(M, x0_3d))),
  t, ts_3d)$

xs_3d : makelist(first(p), p, traj_3d)$
ys_3d : makelist(second(p), p, traj_3d)$
zs_3d : makelist(third(p), p, traj_3d)$

ax_draw3d(
  lines(xs_3d, ys_3d, zs_3d),
  title="3D Spiral Decay",
  xlabel="x", ylabel="y", zlabel="z"
)$
No description has been provided for this image
In [188]:
/* 3D system: time series of each component */
ax_draw2d(
  color="blue", lines(ts_3d, xs_3d), name="x(t)",
  color="red", lines(ts_3d, ys_3d), name="y(t)",
  color="green", lines(ts_3d, zs_3d), name="z(t)",
  title="3D System: Component Time Series",
  xlabel="t", ylabel="value",
  showlegend=true
)$
No description has been provided for this image

Comparing with Symbolic Solution [S+N]¶

For the decoupled 2D system, we can verify that the numeric matrix exponential matches the exact symbolic solution. Since $A$ is diagonal, $e^{At}$ is simply the diagonal matrix of $e^{\lambda_i t}$.

In [189]:
/* Symbolic matrix exponential for the diagonal system */
expAt_sym : matrix([exp(-t), 0], [0, exp(-2*t)])$
print("Symbolic exp(At) =")$
expAt_sym;

/* Verify: x(t) = exp(At) * x0 at t=1 */
x_sym_at_1 : float(subst(t=1, expAt_sym) . matrix([1], [1]));
print("Symbolic x(1) =", x_sym_at_1)$

/* Numeric matrix exponential at t=1 */
M_num : np_expm(np_scale(1.0, A_nd))$
x_num_at_1 : np_to_list(np_matmul(M_num, x0))$
print("Numeric  x(1) =", x_num_at_1)$

print("Agreement confirms np_expm matches the exact solution.")$

Symbolic exp(At) =
matrix([%e^-t,0],[0,%e^-(2*t)])
matrix([0.36787944117144233],[0.1353352832366127])
Symbolic x(1) = matrix([0.36787944117144233],[0.1353352832366127])
Numeric  x(1) = [0.36787944117144233,0.13533528323661267]
Agreement confirms np_expm matches the exact solution.
In [190]:
/* Pointwise comparison across the full trajectory */
sym_xs : makelist(float(exp(-t)), t, ts)$
sym_ys : makelist(float(exp(-2*t)), t, ts)$

/* Compute max absolute error between symbolic and numeric */
err_x : lmax(makelist(abs(xs[i] - sym_xs[i]), i, 1, length(ts)))$
err_y : lmax(makelist(abs(ys[i] - sym_ys[i]), i, 1, length(ts)))$

print("Max error in x(t):", err_x)$
print("Max error in y(t):", err_y)$

Max error in x(t): 1.8041124150158794e-16
Max error in y(t): 2.220446049250313e-16

Multiple Initial Conditions¶

The beauty of the matrix exponential is that we can instantly generate trajectories from any initial condition by computing $e^{At}\mathbf{x}_0$ with different $\mathbf{x}_0$. This reveals the global phase portrait structure.

In [191]:
/* Multiple trajectories for the spiral system from different starting points */
ic_list : [[2.0, 0.0], [0.0, 2.0], [-2.0, 0.0], [0.0, -2.0],
           [1.5, 1.5], [-1.5, 1.5], [-1.5, -1.5], [1.5, -1.5]]$

/* Build a list of line objects, one per initial condition */
plot_cmds : []$
colors : ["#e41a1c", "#377eb8", "#4daf4a", "#984ea3",
          "#ff7f00", "#a65628", "#f781bf", "#999999"]$

for k : 1 thru length(ic_list) do block(
  [ic : ic_list[k], x0_k, traj_k, xk, yk],
  x0_k : ndarray(ic),
  traj_k : makelist(
    block([M : np_expm(np_scale(t, A_spiral_nd))],
      np_to_list(np_matmul(M, x0_k))),
    t, ts_spiral),
  xk : makelist(first(p), p, traj_k),
  yk : makelist(second(p), p, traj_k),
  plot_cmds : append(plot_cmds, [color=colors[k], lines(xk, yk)])
)$

apply(ax_draw2d, append(plot_cmds,
  [title="Damped Spiral: Multiple Initial Conditions",
   xlabel="x", ylabel="y",
   aspect_ratio=true]))$
No description has been provided for this image
In [192]:
/* Same idea for the decoupled system: several initial conditions */
ic_list_2 : [[1.0, 1.0], [-1.0, 1.0], [1.0, -1.0], [-1.0, -1.0],
             [2.0, 0.5], [0.5, 2.0]]$

plot_cmds_2 : []$
colors_2 : ["#e41a1c", "#377eb8", "#4daf4a", "#984ea3",
            "#ff7f00", "#a65628"]$

for k : 1 thru length(ic_list_2) do block(
  [ic : ic_list_2[k], x0_k, traj_k, xk, yk],
  x0_k : ndarray(ic),
  traj_k : makelist(
    block([M : np_expm(np_scale(t, A_nd))],
      np_to_list(np_matmul(M, x0_k))),
    t, ts),
  xk : makelist(first(p), p, traj_k),
  yk : makelist(second(p), p, traj_k),
  plot_cmds_2 : append(plot_cmds_2, [color=colors_2[k], lines(xk, yk)])
)$

apply(ax_draw2d, append(plot_cmds_2,
  [title="Stable Node: Multiple Initial Conditions",
   xlabel="x", ylabel="y",
   aspect_ratio=true]))$
No description has been provided for this image

Summary¶

The matrix exponential $e^{At}$ provides the exact solution to linear ODE systems. This notebook demonstrated:

  1. Numeric computation of $e^{At}$ via np_expm (LAPACK scaling-and-squaring)
  2. Trajectory generation by multiplying $e^{At} \mathbf{x}_0$ at each time step
  3. Verification against symbolic solutions for simple (diagonal) systems
  4. Phase portraits from multiple initial conditions, revealing global dynamics

The eigenvalues of $A$ directly determine the qualitative behaviour:

  • Real negative eigenvalues produce exponential decay (stable node)
  • Complex eigenvalues with negative real part produce spiralling decay
  • The 3D case shows how independent modes combine