EVLRG Example — Companion Matrix Eigenvalues

This example creates the companion matrix for the polynomial \(x^4 - 10x^3 + 35x^2 - 50x + 24\) (roots 1, 2, 3, 4) and computes its eigenvalues using eigensystems.evlrg().

Example Code

"""IMSL EVLRG example: eigenvalues of a companion matrix.

Reproduces the IMSL EVLRG style example:
- Create a 4x4 companion matrix for polynomial x^4 - 10x^3 + 35x^2 - 50x + 24
  whose roots are exactly 1, 2, 3, 4.
- Compute eigenvalues with evlrg().
- Bar chart of real parts of eigenvalues.
- Save to test_output/example_imsl_evlrg.svg

Outputs:
- Printed table of eigenvalues
- SVG bar chart saved to test_output/
- Console output captured to docs/_generated_text/evlrg_output.txt
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from eigensystems import evlrg


def run_demo_evlrg() -> dict:
    """Run IMSL EVLRG example: companion matrix eigenvalues.

    Creates the companion matrix for (x-1)(x-2)(x-3)(x-4) and computes
    its eigenvalues, which should equal the polynomial roots 1, 2, 3, 4.

    Args:
        None

    Returns:
        dict: Result dictionary with keys ``eigenvalues`` (np.ndarray),
            ``sorted_real`` (np.ndarray), and ``plot_path`` (str).
    """
    # Companion matrix for x^4 - 10x^3 + 35x^2 - 50x + 24
    companion = np.array([
        [0.0,  0.0,  0.0, -24.0],
        [1.0,  0.0,  0.0,  50.0],
        [0.0,  1.0,  0.0, -35.0],
        [0.0,  0.0,  1.0,  10.0],
    ])

    result = evlrg(companion)
    sorted_real = np.sort(result.eigenvalues.real)

    print("\nIMSL EVLRG Example: Companion Matrix Eigenvalues")
    print("Polynomial: x^4 - 10x^3 + 35x^2 - 50x + 24")
    print("Expected roots: 1, 2, 3, 4")
    print("-" * 45)
    print(f"{'Index':<8} {'Eigenvalue (complex)':>25}")
    print("-" * 45)
    for i, ev in enumerate(sorted(_sorted_complex(result.eigenvalues),
                                   key=lambda z: z.real)):
        print(f"  {i+1:<6} {ev.real:>12.6f} + {ev.imag:+.6f}i")
    print("-" * 45)
    print(f"Sorted real parts: {sorted_real}")

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

    fig, ax = plt.subplots(figsize=(7, 5))
    idx = np.arange(1, len(sorted_real) + 1)
    ax.bar(idx, sorted_real, color="#be123c", alpha=0.85, edgecolor="#7f1d1d")
    ax.set_xticks(idx)
    ax.set_xticklabels([f"λ{i}" for i in idx])
    ax.set_xlabel("Eigenvalue index")
    ax.set_ylabel("Real part of eigenvalue")
    ax.set_title("IMSL EVLRG: Companion Matrix Eigenvalues\n"
                 r"$x^4 - 10x^3 + 35x^2 - 50x + 24$, roots = {1, 2, 3, 4}")
    ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved: {plot_path}")

    return {
        "eigenvalues": result.eigenvalues,
        "sorted_real": sorted_real,
        "plot_path": str(plot_path),
    }


def _sorted_complex(arr: np.ndarray) -> np.ndarray:
    """Sort a complex array by real part ascending.

    Args:
        arr (np.ndarray): Array of complex numbers.

    Returns:
        np.ndarray: Sorted array.
    """
    return arr[np.argsort(arr.real)]


if __name__ == "__main__":
    run_demo_evlrg()

Plot Output

Real parts of the four eigenvalues. The companion matrix eigenvalues equal the polynomial roots: 1, 2, 3, 4.

Console Output

IMSL EVLRG Example: Companion Matrix Eigenvalues
Polynomial: x^4 - 10x^3 + 35x^2 - 50x + 24
Expected roots: 1, 2, 3, 4
---------------------------------------------
Index         Eigenvalue (complex)
---------------------------------------------
  1          1.000000 + +0.000000i
  2          2.000000 + +0.000000i
  3          3.000000 + +0.000000i
  4          4.000000 + +0.000000i
---------------------------------------------
Sorted real parts: [1. 2. 3. 4.]

Plot saved: test_output\example_imsl_evlrg.svg