ZPLRC Example — Polynomial Root-Finding

This example finds all roots of the real polynomial \(p(x) = x^4 - 5x^3 + 5x^2 + 5x - 6 = (x+1)(x-1)(x-2)(x-3)\) using the eigenvalue method (IMSL ZPLRC). The roots are located in the complex plane and compared against the known exact values \(\{-1, 1, 2, 3\}\).

Example Code

"""IMSL ZPLRC example: roots of a real polynomial using eigenvalue method.

Reproduces the IMSL ZPLRC published example:
- Polynomial: p(x) = x^4 - 5x^3 + 5x^2 + 5x - 6 = (x-1)(x-2)(x+1)(x-3)
- Expected roots: approximately {-1, 1, 2, 3}

Outputs:
- Table printed to stdout
- SVG scatter plot of roots in the complex plane saved to test_output/demo_imsl_zplrc.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from nonlinearequations import polynomial_roots_real


def run_demo_imsl_zplrc() -> Dict[str, object]:
    """Run IMSL ZPLRC example: roots of (x+1)(x-1)(x-2)(x-3).

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``roots`` (np.ndarray),
            ``fval`` (float), and ``plot_path`` (str).
    """
    # coefficients ascending: p(x) = -6 + 5x + 5x^2 - 5x^3 + x^4
    coeffs = np.array([-6.0, 5.0, 5.0, -5.0, 1.0])
    result = polynomial_roots_real(coeffs)
    roots = result.x
    roots_sorted = roots[np.argsort(roots.real)]

    print("\nIMSL ZPLRC Example: Roots of x^4 - 5x^3 + 5x^2 + 5x - 6")
    print("-" * 65)
    print(f"{'#':<4} {'Real Part':>16} {'Imag Part':>16} {'|p(root)|':>14}")
    print("-" * 65)
    poly_desc = coeffs[::-1]
    for i, z in enumerate(roots_sorted):
        resid = abs(np.polyval(poly_desc, z))
        print(f"{i+1:<4} {z.real:>16.10f} {z.imag:>16.2e} {resid:>14.2e}")
    print("-" * 65)
    print(f"{'Max residual':<30} {result.fval:>14.2e}")
    print(f"{'Converged':<30} {str(result.success):>14}")
    print(f"{'Expected roots':<30} {'[-1, 1, 2, 3]':>14}")
    print("-" * 65)

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

    # Plot roots in complex plane
    theta = np.linspace(0, 2 * np.pi, 300)
    unit_x = np.cos(theta)
    unit_y = np.sin(theta)

    fig, ax = plt.subplots(figsize=(7, 6))
    ax.plot(unit_x, unit_y, color="#aaaaaa", linewidth=1.0, linestyle="--", label="Unit circle")
    ax.scatter(
        roots.real, roots.imag,
        color="#0e7490", s=90, zorder=5, label="Computed roots"
    )
    expected = np.array([-1.0, 1.0, 2.0, 3.0])
    ax.scatter(expected, np.zeros_like(expected), color="#d62728", s=40, marker="x", zorder=6, label="Expected roots")
    for z in roots_sorted:
        ax.annotate(
            f"{z.real:.4f}", xy=(z.real, z.imag),
            xytext=(z.real, z.imag + 0.15),
            ha="center", fontsize=9, color="#0e7490"
        )
    ax.axhline(0, color="#888888", linewidth=0.7)
    ax.axvline(0, color="#888888", linewidth=0.7)
    ax.set_xlabel("Real")
    ax.set_ylabel("Imaginary")
    ax.set_title("IMSL ZPLRC: Roots of $x^4 - 5x^3 + 5x^2 + 5x - 6$")
    ax.legend()
    ax.grid(True, alpha=0.25)
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"roots": roots, "fval": result.fval, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_zplrc()

Plot Output

Roots of \(x^4 - 5x^3 + 5x^2 + 5x - 6\) plotted in the complex plane. Computed roots (blue circles) coincide with the expected values \(-1, 1, 2, 3\) (red crosses) on the real axis.

Console Output

IMSL ZPLRC Example: Roots of x^4 - 5x^3 + 5x^2 + 5x - 6
-----------------------------------------------------------------
#           Real Part        Imag Part      |p(root)|
-----------------------------------------------------------------
1       -1.0000000000         0.00e+00       5.33e-15
2        1.0000000000         0.00e+00       4.44e-15
3        2.0000000000         0.00e+00       1.87e-14
4        3.0000000000        -0.00e+00       4.97e-14
-----------------------------------------------------------------
Max residual                         4.97e-14
Converged                                True
Expected roots                  [-1, 1, 2, 3]
-----------------------------------------------------------------