Stability Analysis and BVP Shooting¶
This example demonstrates three numerical analysis utilities:
stability_check()— eigenvalues of the numerical Jacobian at equilibrium points.ivp_jacobian()— numerical Jacobian of \(f\) evaluated along a trajectory.bvpms()— single-shooting boundary value problem solver.
Part A — Pendulum Equilibria¶
The nonlinear pendulum \(\theta' = \omega\), \(\omega' = -\sin\theta\) has two equilibria in \([0, 2\pi]\):
\((\theta, \omega) = (0, 0)\): stable centre with pure-imaginary eigenvalues \(\pm i\).
\((\theta, \omega) = (\pi, 0)\): saddle point with real eigenvalues \(\pm 1\) (unstable).
Part B — Jacobian Along a Trajectory¶
ivp_jacobian evaluates \(\partial f / \partial y\) at eight equally
spaced points on a small-amplitude trajectory, tracking how the linearised
dynamics vary with the angle \(\theta(t)\).
Part C — BVP Shooting¶
The BVP
\[y'' + y = 0, \quad y(0) = 0, \quad y(\pi/2) = 1\]
has the exact solution \(y(x) = \sin x\). bvpms uses Brent’s root
finding algorithm to locate the initial slope \(y'(0)\) that satisfies
the right boundary condition, integrating the IVP with ivprk.
Example Code¶
"""IMSL stability analysis and BVP shooting example.
Demonstrates ``stability_check``, ``ivp_jacobian``, and ``bvpms``:
Part A — Pendulum equilibria:
θ' = ω, ω' = -sin(θ)
Stable fixed point at (0, 0), unstable at (π, 0).
``stability_check`` is called at both equilibria.
Part B — Numerical Jacobian along a trajectory:
``ivp_jacobian`` is evaluated at several points on a small-amplitude
pendulum trajectory to track how the local linearisation evolves.
Part C — BVP shooting method:
y'' + y = 0, y(0) = 0, y(π/2) = 1
Exact solution: y(x) = sin(x).
``bvpms`` locates the initial slope that satisfies the right boundary
condition.
Outputs:
- Eigenvalue summary printed to stdout.
- SVG plot saved to ``test_output/example_imsl_stability.svg``.
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from differentialequations import bvpms, ivp_jacobian, ivprk, stability_check
# ---------------------------------------------------------------------------
# Part A: Pendulum equilibria
# ---------------------------------------------------------------------------
def pendulum(t: float, y: np.ndarray) -> list[float]:
"""Nonlinear pendulum ODE right-hand side.
Args:
t (float): Time (unused — autonomous system).
y (np.ndarray): State ``[θ, ω]``.
Returns:
list[float]: Derivatives ``[ω, -sin(θ)]``.
"""
theta, omega = y
return [omega, -np.sin(theta)]
def run_demo_imsl_stability() -> Dict[str, object]:
"""Run stability analysis, Jacobian computation, and BVP shooting examples.
Part A computes eigenvalues of the Jacobian at the pendulum equilibria
(stable and unstable) using ``stability_check``.
Part B evaluates ``ivp_jacobian`` at several points along a trajectory of
the small-amplitude pendulum to illustrate how the Jacobian varies.
Part C solves the BVP y'' + y = 0, y(0)=0, y(π/2)=1 using ``bvpms``.
Args:
None
Returns:
Dict[str, object]: Keys ``eigenvalues_stable``, ``eigenvalues_unstable``
(complex ndarrays), ``bvp_t``, ``bvp_y`` (ndarrays), and
``plot_path`` (str).
"""
# -----------------------------------------------------------------------
# Part A: stability_check at stable (0,0) and unstable (π,0) equilibria
# -----------------------------------------------------------------------
y_stable = [0.0, 0.0]
y_unstable = [np.pi, 0.0]
eigs_stable = stability_check(pendulum, y_stable, t=0.0)
eigs_unstable = stability_check(pendulum, y_unstable, t=0.0)
print("\nIMSL Stability Analysis: Nonlinear Pendulum")
print("=" * 60)
print("\nPart A — Equilibrium eigenvalues (stability_check)")
print(f" Stable fixed point (θ=0, ω=0): eigenvalues = {eigs_stable}")
print(f" Re parts: {eigs_stable.real}")
print(f" → All Re<0: {np.all(eigs_stable.real < 0)}")
print(f"\n Unstable fixed point (θ=π, ω=0): eigenvalues = {eigs_unstable}")
print(f" Re parts: {eigs_unstable.real}")
print(f" → Any Re>0 (unstable): {np.any(eigs_unstable.real > 0)}")
# -----------------------------------------------------------------------
# Part B: ivp_jacobian along a small-amplitude trajectory
# -----------------------------------------------------------------------
print("\nPart B — Jacobian along trajectory (ivp_jacobian)")
print(f" {'t':>6} {'J[0,0]':>10} {'J[0,1]':>10} {'J[1,0]':>10} {'J[1,1]':>10}")
print(" " + "-" * 56)
traj = ivprk(pendulum, (0.0, 2.0 * np.pi), [0.3, 0.0])
sample_indices = np.linspace(0, traj.t.size - 1, 8, dtype=int)
jac_samples = []
for idx in sample_indices:
ti = traj.t[idx]
yi = traj.y[:, idx]
J = ivp_jacobian(pendulum, ti, yi)
jac_samples.append((ti, J))
print(
f" {ti:>6.3f} {J[0, 0]:>10.4f} {J[0, 1]:>10.4f} "
f"{J[1, 0]:>10.4f} {J[1, 1]:>10.4f}"
)
# -----------------------------------------------------------------------
# Part C: BVP shooting — y'' + y = 0, y(0)=0, y(π/2)=1
# -----------------------------------------------------------------------
print("\nPart C — BVP shooting method (bvpms)")
print(" y'' + y = 0, y(0) = 0, y(π/2) = 1")
print(" Exact solution: y(x) = sin(x)")
def bvp_ode(x: float, y: np.ndarray) -> list[float]:
"""Second-order ODE y''+y=0 as first-order system.
Args:
x (float): Independent variable.
y (np.ndarray): State ``[y, y']``.
Returns:
list[float]: Derivatives ``[y', -y]``.
"""
return [y[1], -y[0]]
def bvp_bc(ya: np.ndarray, yb: np.ndarray) -> np.ndarray:
"""Boundary conditions y(0)=0, y(π/2)=1.
Args:
ya (np.ndarray): State at left boundary.
yb (np.ndarray): State at right boundary.
Returns:
np.ndarray: Boundary residuals.
"""
return np.array([ya[0], yb[0] - 1.0])
x_mesh = np.linspace(0.0, np.pi / 2.0, 30)
# Initial guess: linear ramp consistent with left BC
y_guess = np.zeros((2, x_mesh.size))
y_guess[0] = x_mesh / (np.pi / 2.0)
y_guess[1] = 1.0
bvp_res = bvpms(bvp_ode, bvp_bc, x_mesh, y_guess, tol=1e-4)
y_exact = np.sin(bvp_res.t)
max_err = np.max(np.abs(bvp_res.y[0] - y_exact))
print(f"\n {'x':>8} {'y_shoot':>12} {'y_exact':>12} {'|error|':>12}")
print(" " + "-" * 52)
display_idx = np.linspace(0, bvp_res.t.size - 1, 8, dtype=int)
for idx in display_idx:
xi = bvp_res.t[idx]
ys = bvp_res.y[0, idx]
ye = np.sin(xi)
print(f" {xi:>8.4f} {ys:>12.6f} {ye:>12.6f} {abs(ys - ye):>12.2e}")
print(f"\n Max error vs sin(x): {max_err:.2e}")
print(f" Solver message: {bvp_res.message}")
print("=" * 60)
# -----------------------------------------------------------------------
# Plots
# -----------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_stability.svg"
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Panel 1: eigenvalues of Jacobian on complex plane
ax = axes[0]
ax.axhline(0, color="k", linewidth=0.8)
ax.axvline(0, color="k", linewidth=0.8)
ax.scatter(
eigs_stable.real,
eigs_stable.imag,
s=120,
marker="o",
color="#0369a1",
zorder=5,
label="Stable (θ=0): Re<0",
)
ax.scatter(
eigs_unstable.real,
eigs_unstable.imag,
s=120,
marker="^",
color="#dc2626",
zorder=5,
label="Unstable (θ=π): Re>0",
)
for ev in eigs_stable:
ax.annotate(
f"({ev.real:.2f}, {ev.imag:.2f}i)",
(ev.real, ev.imag),
textcoords="offset points",
xytext=(6, 6),
fontsize=7,
color="#0369a1",
)
for ev in eigs_unstable:
ax.annotate(
f"({ev.real:.2f}, {ev.imag:.2f}i)",
(ev.real, ev.imag),
textcoords="offset points",
xytext=(6, -12),
fontsize=7,
color="#dc2626",
)
ax.set_xlabel("Re(λ)")
ax.set_ylabel("Im(λ)")
ax.set_title("Eigenvalues at Equilibria\n(stability_check)")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# Panel 2: J[1,0] = -cos(θ) along trajectory (dominant varying entry)
ax = axes[1]
t_vals = [s[0] for s in jac_samples]
j10_vals = [s[1][1, 0] for s in jac_samples]
ax.plot(traj.t, traj.y[0], color="#0369a1", linewidth=1.5, label="θ(t)")
ax2 = ax.twinx()
ax2.plot(t_vals, j10_vals, "ro--", linewidth=1.5, markersize=6, label="J[1,0]")
ax2.set_ylabel("J[1,0] = −cos(θ)", color="#dc2626")
ax2.tick_params(axis="y", labelcolor="#dc2626")
ax.set_xlabel("t")
ax.set_ylabel("θ (rad)", color="#0369a1")
ax.tick_params(axis="y", labelcolor="#0369a1")
ax.set_title("Jacobian Along Trajectory\n(ivp_jacobian)")
lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, fontsize=8, loc="upper right")
ax.grid(True, alpha=0.3)
# Panel 3: BVP shooting vs exact solution
ax = axes[2]
x_fine = np.linspace(0.0, np.pi / 2.0, 200)
ax.plot(x_fine, np.sin(x_fine), "k--", linewidth=1.5, label="Exact: sin(x)")
ax.plot(
bvp_res.t,
bvp_res.y[0],
"o-",
color="#16a34a",
linewidth=2.0,
markersize=4,
label=f"bvpms (max err={max_err:.1e})",
)
ax.set_xlabel("x")
ax.set_ylabel("y(x)")
ax.set_title("BVP Shooting Method (bvpms)\ny''+y=0, y(0)=0, y(π/2)=1")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.suptitle(
"IMSL Stability & BVP: stability_check, ivp_jacobian, bvpms",
fontsize=12,
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"eigenvalues_stable": eigs_stable,
"eigenvalues_unstable": eigs_unstable,
"bvp_t": bvp_res.t,
"bvp_y": bvp_res.y,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_stability()
Input (Console)¶
Run the script from the package root:
python examples/example_imsl_stability.py
Console Output¶
IMSL Stability Analysis: Nonlinear Pendulum
============================================================
Part A ΓÇö Equilibrium eigenvalues (stability_check)
Stable fixed point (θ=0, ω=0): eigenvalues = [0.+1.j 0.-1.j]
Re parts: [0. 0.]
→ All Re<0: False
Unstable fixed point (θ=π, ω=0): eigenvalues = [ 1. -1.]
Re parts: [ 1. -1.]
→ Any Re>0 (unstable): True
Part B ΓÇö Jacobian along trajectory (ivp_jacobian)
t J[0,0] J[0,1] J[1,0] J[1,1]
--------------------------------------------------------
0.000 0.0000 1.0000 -0.9553 0.0000
0.354 0.0000 1.0000 -0.9606 0.0000
1.478 0.0000 1.0000 -0.9995 0.0000
2.428 0.0000 1.0000 -0.9750 0.0000
3.287 0.0000 1.0000 -0.9560 0.0000
4.366 0.0000 1.0000 -0.9941 0.0000
5.320 0.0000 1.0000 -0.9866 0.0000
6.283 0.0000 1.0000 -0.9554 0.0000
Part C ΓÇö BVP shooting method (bvpms)
y'' + y = 0, y(0) = 0, y(π/2) = 1
Exact solution: y(x) = sin(x)
x y_shoot y_exact |error|
----------------------------------------------------
0.0000 0.000000 0.000000 0.00e+00
0.2167 0.214971 0.214970 2.49e-07
0.4333 0.419885 0.419889 3.66e-06
0.6500 0.605165 0.605174 9.48e-06
0.8666 0.762161 0.762162 1.20e-06
1.0833 0.883516 0.883512 3.52e-06
1.3000 0.963507 0.963550 4.26e-05
1.5708 1.000001 1.000000 8.79e-07
Max error vs sin(x): 5.28e-05
Solver message: Shooting method converged
============================================================
Plot saved to: test_output\example_imsl_stability.svg
Plot Output¶
Generated SVG plot:
Left: eigenvalues of the Jacobian at both equilibria plotted on the complex
plane — the stable centre has pure-imaginary eigenvalues (on the imaginary
axis) while the saddle has real eigenvalues of opposite sign. Centre:
the \(J_{10}\) entry (\(-\cos\theta\)) versus time, overlaid on
the trajectory. Right: bvpms shooting solution compared against the
exact \(\sin x\).¶
Description¶
stability_check internally calls ivp_jacobian to build the finite-
difference Jacobian and then returns its eigenvalues. The pendulum at
\(\theta = 0\) is a Lyapunov-stable centre (pure-imaginary eigenvalues,
conservative system), while at \(\theta = \pi\) it is an unstable saddle
(one positive real eigenvalue). The BVP shooting method reaches a maximum
error of order \(10^{-5}\) against the exact solution with default
tolerances.