Symbolic Jacobian and Stability Analysis¶

This notebook demonstrates the symbolic-numeric bridge for dynamical systems analysis. We use Maxima's diff() to derive Jacobian matrices exactly, solve() to find equilibrium points, and np_eig to compute eigenvalues numerically for stability classification.

The workflow is: symbolic differentiation -> symbolic algebra -> numeric eigenvalues -> visualization. Each step leverages the right tool for the job.

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

The Jacobian [S+N]¶

Consider a predator-prey variant:

$$\dot{x} = x(1-x) - xy \qquad \dot{y} = y(x - 1)$$

Here $x$ is the prey population (with logistic self-limitation) and $y$ is the predator. We compute the Jacobian matrix symbolically using diff().

In [98]:
/* Define the system */
f1 : x*(1 - x) - x*y$
f2 : y*(x - 1)$

/* Compute the Jacobian symbolically */
J : matrix(
  [diff(f1, x), diff(f1, y)],
  [diff(f2, x), diff(f2, y)]
);

print("Jacobian:")$
J;

matrix([-y-2*x+1,-x],[y,x-1])
Jacobian:
Out[98]:

Fixed Points [S+N]¶

Equilibria are where both $\dot{x} = 0$ and $\dot{y} = 0$ simultaneously. Maxima's solve() finds them exactly.

In [99]:
/* Find all equilibrium points */
fps : solve([f1 = 0, f2 = 0], [x, y]);

print("Fixed points:")$
for fp in fps do print(fp)$

[[x = 0,y = 0],[x = 1,y = 0]]
Fixed points:
[x = 0,y = 0]
[x = 1,y = 0]

Stability Classification [S+N]¶

For each fixed point, we substitute into the Jacobian and compute eigenvalues. The eigenvalues determine stability:

  • Both eigenvalues have negative real part: stable (node or spiral)
  • Both positive: unstable
  • Mixed signs: saddle point

We use Maxima's eigenvalues() for the symbolic case and np_eig for the numeric case.

In [100]:
/* Evaluate Jacobian at each fixed point and classify */
for fp in fps do block(
  [J_fp, ev_list, re1, re2, cl],
  J_fp : subst(fp, J),
  ev_list : first(eigenvalues(J_fp)),
  re1 : float(realpart(ev_list[1])),
  re2 : float(realpart(ev_list[2])),
  cl : if is(re1 < 0 and re2 < 0) then "STABLE"
       elseif is(re1 > 0 and re2 > 0) then "UNSTABLE"
       elseif is(re1 * re2 < 0) then "SADDLE"
       else "MARGINAL / CENTRE",
  print("--- Fixed point:", fp, "---"),
  print("  J =", J_fp),
  print("  Eigenvalues =", ev_list),
  print("  Classification:", cl)
)$

--- Fixed point: [x = 0,y = 0] ---
  J = matrix([1,0],[0,-1])
  Eigenvalues = [-1,1]
  Classification: SADDLE
--- Fixed point: [x = 1,y = 0] ---
  J = matrix([-1,-1],[0,0])
  Eigenvalues = [-1,0]
  Classification: MARGINAL / CENTRE

Numeric Eigenvalues with np_eig¶

For more complex systems where symbolic eigenvalues may not simplify, we can evaluate the Jacobian numerically and use np_eig from the numerics package (wrapping LAPACK's dgeev).

In [101]:
/* Numeric eigenvalue computation at (1, 0) using np_eig */
J_at_10 : float(subst([x=1, y=0], J))$
print("J at (1,0) =", J_at_10)$

A_np : ndarray(J_at_10)$
[vals, vecs] : np_eig(A_np)$
print("Eigenvalues (numeric):", np_to_list(vals))$
print("Eigenvectors:")$
np_to_matrix(vecs);

J at (1,0) = matrix([-1.0,-1.0],[0.0,0.0])
Eigenvalues (numeric): [-1.0,0.0]
Eigenvectors:
Out[101]:

Phase Portrait¶

The vector field shows the flow of the dynamical system at every point. Streamlines trace trajectories from various initial conditions, revealing the global behaviour around the fixed points.

In [102]:
/* Vector field of the predator-prey system */
ax_draw2d(
  color="#cccccc", ngrid=12,
  ax_vector_field(f1, f2, x, -0.5, 2.5, y, -0.5, 2.5),
  title="Predator-Prey: Vector Field",
  xlabel="x (prey)", ylabel="y (predator)",
  aspect_ratio=true
)$
No description has been provided for this image
In [103]:
/* Phase portrait: vector field with streamlines */
ax_draw2d(
  color="#dddddd", ngrid=12,
  ax_vector_field(f1, f2, x, -0.5, 2.5, y, -0.5, 2.5),
  color="#e41a1c",
  initial_points=[[0.1, 0.1], [0.5, 0.5], [1.5, 0.2],
                  [0.2, 1.0], [2.0, 0.5], [0.8, 1.5]],
  ax_streamline(f1, f2, x, -0.5, 2.5, y, -0.5, 2.5),
  title="Predator-Prey: Phase Portrait",
  xlabel="x (prey)", ylabel="y (predator)",
  aspect_ratio=true
)$
No description has been provided for this image

Van der Pol Oscillator [S+N]¶

The Van der Pol oscillator is a classic nonlinear system with a limit cycle:

$$\dot{x} = y \qquad \dot{y} = \mu(1 - x^2)y - x$$

For $\mu > 0$ the origin is an unstable spiral, and trajectories converge to a stable limit cycle. We derive the Jacobian symbolically, compute eigenvalues at the origin, and visualize the limit cycle.

In [104]:
/* Van der Pol system with mu = 1 */
mu : 1$
g1 : y$
g2 : mu*(1 - x^2)*y - x$

/* Symbolic Jacobian */
J_vdp : matrix(
  [diff(g1, x), diff(g1, y)],
  [diff(g2, x), diff(g2, y)]
);

print("Van der Pol Jacobian:")$
J_vdp;

matrix([0,1],[-(2*x*y)-1,1-x^2])
Van der Pol Jacobian:
Out[104]:
In [105]:
/* Evaluate at the origin — the only fixed point */
J_vdp_origin : subst([x=0, y=0], J_vdp)$
print("J at origin:", J_vdp_origin)$

/* Symbolic eigenvalues */
evals_vdp : eigenvalues(J_vdp_origin);
print("Eigenvalues:", evals_vdp[1])$

/* Numeric check via np_eig */
A_vdp : ndarray(float(J_vdp_origin))$
[vals_vdp, vecs_vdp] : np_eig(A_vdp)$
print("Eigenvalues (numeric):", np_to_list(vals_vdp))$

J at origin: matrix([0,1],[-1,1])
[[-((sqrt(3)*%i-1)/2),(sqrt(3)*%i+1)/2],[1,1]]
Eigenvalues: [-((sqrt(3)*%i-1)/2),(sqrt(3)*%i+1)/2]
Eigenvalues (numeric):
                      [0.8660254037844386*%i+0.49999999999999994,
                       0.49999999999999994-0.8660254037844386*%i]

The eigenvalues at the origin have positive real part ($\mu/2 > 0$), confirming the origin is an unstable spiral. All nearby trajectories spiral outward — but they are captured by the stable limit cycle.

In [106]:
/* Phase portrait of the Van der Pol oscillator */
ax_draw2d(
  color="#dddddd", ngrid=12,
  ax_vector_field(g1, g2, x, -4, 4, y, -4, 4),
  color="#377eb8",
  initial_points=[[0.1, 0.0], [3.5, 0.0], [0.0, 0.5],
                  [-2.0, 2.0], [1.0, -3.0]],
  t_range=[0, 30],
  ax_streamline(g1, g2, x, -4, 4, y, -4, 4),
  title="Van der Pol Oscillator (mu=1): Limit Cycle",
  xlabel="x", ylabel="y",
  aspect_ratio=true
)$
No description has been provided for this image

Varying the Damping Parameter¶

As $\mu$ increases, the limit cycle becomes more pronounced and the oscillation becomes more "relaxation-like" (sharp transitions). We compare the Jacobian eigenvalues at the origin for different $\mu$.

In [107]:
/* Eigenvalues at origin for varying mu — symbolic formula */
J_vdp_general : matrix(
  [0, 1],
  [-1, mu_var]
)$
evals_general : eigenvalues(J_vdp_general);

print("General eigenvalues at origin:")$
print(evals_general[1])$

/* Evaluate for specific mu values */
for m in [0.1, 0.5, 1.0, 2.0, 5.0] do block(
  [ev : float(subst(mu_var=m, evals_general[1]))],
  print(sconcat("mu = ", m, ":"), ev)
)$

[[-((sqrt(mu_var^2-4)-mu_var)/2),(sqrt(mu_var^2-4)+mu_var)/2],[1,1]]
General eigenvalues at origin:
[-((sqrt(mu_var^2-4)-mu_var)/2),(sqrt(mu_var^2-4)+mu_var)/2]
mu = 0.1: [-(0.5*(1.997498435543818*%i-0.1)),0.5*(1.997498435543818*%i+0.1)]
mu = 0.5: [-(0.5*(1.9364916731037085*%i-0.5)),0.5*(1.9364916731037085*%i+0.5)]
mu = 1.0: [-(0.5*(1.7320508075688772*%i-1.0)),0.5*(1.7320508075688772*%i+1.0)]
mu = 2.0: [1.0,1.0]
mu = 5.0: [0.20871215252208009,4.7912878474779195]
In [108]:
/* Phase portrait with mu = 3: strongly nonlinear limit cycle */
g2_strong : 3*(1 - x^2)*y - x$

ax_draw2d(
  color="#dddddd", ngrid=12,
  ax_vector_field(g1, g2_strong, x, -4, 4, y, -6, 6),
  color="#e41a1c",
  initial_points=[[0.1, 0.0], [3.5, 0.0], [-2.0, 3.0]],
  t_range=[0, 40],
  ax_streamline(g1, g2_strong, x, -4, 4, y, -6, 6),
  title="Van der Pol (mu=3): Relaxation Oscillation",
  xlabel="x", ylabel="y"
)$
No description has been provided for this image

Summary¶

Key insight: Maxima's symbolic differentiation gives us the exact Jacobian matrix for any smooth system. No finite-difference approximation, no numerical error in the derivatives. We then:

  1. Solve for equilibria symbolically with solve()
  2. Substitute fixed-point coordinates into the Jacobian
  3. Compute eigenvalues — symbolically with eigenvalues() for simple cases, or numerically with np_eig (LAPACK) for larger systems
  4. Classify stability: stable / unstable / saddle / spiral
  5. Visualise the phase portrait with ax_vector_field and ax_streamline

This is the symbolic-numeric bridge at work: exact analysis where it matters (derivatives, algebra), numerical muscle where it's needed (eigenvalues, trajectory integration, plotting).