Phase Portrait Example — Pendulum & Lotka-Volterra

This example uses differentialequations.phase_portrait() to generate multi-trajectory phase portraits for two classical nonlinear ODE systems.

a) Simple pendulum

\[\begin{split}\theta'' = -\sin\theta \quad\Longrightarrow\quad \begin{bmatrix}\theta \\ \omega\end{bmatrix}' = \begin{bmatrix}\omega \\ -\sin\theta\end{bmatrix}\end{split}\]

b) Lotka-Volterra predator-prey (\(a=1.5,\;b=1,\;c=3,\;d=1\))

\[x' = a\,x - b\,x\,y, \quad y' = d\,x\,y - c\,y\]

Example Code

"""IMSL phase portrait example: pendulum and Lotka-Volterra.

Demonstrates the use of phase_portrait() to generate multi-trajectory phase
portraits for two classical ODE systems:

    a) Simple pendulum:  θ'' = -sin(θ)  →  [θ, ω]' = [ω, -sin(θ)]
    b) Lotka-Volterra predator-prey:
           x' = a*x - b*x*y   (prey,     a=1.5, b=1)
           y' = d*x*y - c*y   (predator, d=1,   c=3)

Outputs:
- Summary printed to stdout
- Two-panel SVG saved to test_output/example_imsl_phase_portrait.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 phase_portrait


def _pendulum_rhs(t: float, y: np.ndarray) -> list:
    """Simple pendulum right-hand side.

    Args:
        t (float): Time (unused, autonomous).
        y (np.ndarray): State [θ, ω].

    Returns:
        list: [dθ/dt, dω/dt] = [ω, -sin(θ)].
    """
    return [y[1], -np.sin(y[0])]


def _lotka_volterra_rhs(t: float, y: np.ndarray) -> list:
    """Lotka-Volterra predator-prey right-hand side.

    Args:
        t (float): Time (unused, autonomous).
        y (np.ndarray): State [x, y] = [prey, predator].

    Returns:
        list: [dx/dt, dy/dt].
    """
    a, b, c, d = 1.5, 1.0, 3.0, 1.0
    x, yp = y
    return [a * x - b * x * yp, d * x * yp - c * yp]


def run_demo_imsl_phase_portrait() -> Dict[str, object]:
    """Run IMSL phase portrait example for pendulum and Lotka-Volterra.

    Uses phase_portrait() with multiple initial conditions to generate phase
    trajectories for both systems.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``pendulum_results``,
            ``lotka_results``, and ``plot_path`` (str).
    """
    # --- Pendulum initial conditions ---
    angles = np.linspace(0.2, 2.8, 6)
    pendulum_ics = [[th, 0.0] for th in angles]
    pendulum_ics += [[-2.5, 1.5], [2.5, -1.5]]
    t_eval_pend = np.linspace(0.0, 15.0, 600)
    pendulum_results = phase_portrait(
        _pendulum_rhs, (0.0, 15.0), pendulum_ics, t_eval=t_eval_pend
    )

    # --- Lotka-Volterra initial conditions ---
    prey_vals = [0.5, 1.0, 1.5, 2.0, 2.5]
    lv_ics = [[x0, 1.0] for x0 in prey_vals]
    t_eval_lv = np.linspace(0.0, 12.0, 800)
    lotka_results = phase_portrait(
        _lotka_volterra_rhs, (0.0, 12.0), lv_ics, t_eval=t_eval_lv
    )

    print("\nIMSL Phase Portrait Example")
    print("=" * 55)
    print(f"  Pendulum ICs    : {len(pendulum_ics)} trajectories")
    print(f"  Lotka-Volterra  : {len(lv_ics)} trajectories  (a=1.5, b=1, c=3, d=1)")
    print()
    print("  Pendulum equilibria: (0, 0) — stable centre")
    print("                       (±π, 0) — saddle (separatrix)")
    print()
    print("  Lotka-Volterra interior equilibrium: (c/d, a/b) = (3.0, 1.5)")
    print("=" * 55)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_phase_portrait.svg"

    colors_pend = plt.cm.cool(np.linspace(0.1, 0.9, len(pendulum_results)))
    colors_lv = plt.cm.viridis(np.linspace(0.1, 0.9, len(lotka_results)))

    fig, axes = plt.subplots(1, 2, figsize=(13, 5))

    # Pendulum phase portrait
    for res, col in zip(pendulum_results, colors_pend):
        axes[0].plot(res.y[0], res.y[1], color=col, linewidth=1.5, alpha=0.85)
    axes[0].axhline(0, color="black", linewidth=0.6, linestyle=":")
    axes[0].axvline(0, color="black", linewidth=0.6, linestyle=":")
    axes[0].scatter([0, np.pi, -np.pi], [0, 0, 0], s=40, color="#e11d48",
                    zorder=5, label="Equilibria")
    axes[0].set_xlabel(r"$\theta$ (rad)")
    axes[0].set_ylabel(r"$\omega$ (rad/s)")
    axes[0].set_title("Simple Pendulum Phase Portrait")
    axes[0].legend(fontsize=9)
    axes[0].grid(True, alpha=0.25)

    # Lotka-Volterra phase portrait
    eq_x, eq_y = 3.0, 1.5
    for res, col in zip(lotka_results, colors_lv):
        axes[1].plot(res.y[0], res.y[1], color=col, linewidth=1.5, alpha=0.85)
    axes[1].scatter([eq_x], [eq_y], s=60, color="#e11d48", zorder=5,
                    label=f"Equil. ({eq_x}, {eq_y})")
    axes[1].set_xlabel("Prey population  x")
    axes[1].set_ylabel("Predator population  y")
    axes[1].set_title("Lotka-Volterra Phase Portrait")
    axes[1].legend(fontsize=9)
    axes[1].grid(True, alpha=0.25)

    fig.suptitle("IMSL phase_portrait: Pendulum & Lotka-Volterra", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {
        "pendulum_results": pendulum_results,
        "lotka_results": lotka_results,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_phase_portrait()

Plot Output

Left: pendulum phase portrait with closed orbits around the stable centre at the origin and separatrices through \((\pm\pi, 0)\). Right: Lotka-Volterra closed orbits around the interior equilibrium \((c/d,\, a/b) = (3.0,\, 1.5)\).

Console Output

IMSL Phase Portrait Example
=======================================================
  Pendulum ICs    : 8 trajectories
  Lotka-Volterra  : 5 trajectories  (a=1.5, b=1, c=3, d=1)

  Pendulum equilibria: (0, 0) ΓÇö stable centre
                       (±π, 0) — saddle (separatrix)

  Lotka-Volterra interior equilibrium: (c/d, a/b) = (3.0, 1.5)
=======================================================