NEQNF Example — Nonlinear System Solver

This example solves a 2-variable nonlinear system using a finite-difference Jacobian (IMSL NEQNF):

\[ \begin{align}\begin{aligned}x_1^2 + x_2^2 = 1 \quad \text{(unit circle)}\\x_1 - x_2 = 0 \quad \text{(diagonal line)}\end{aligned}\end{align} \]

The intersection gives two solutions: \((1/\sqrt{2},\, 1/\sqrt{2})\) and \((-1/\sqrt{2},\, -1/\sqrt{2})\). Starting from \((0.8, 0.6)\), the solver converges to the first solution.

Example Code

"""IMSL NEQNF example: solve a 2-variable nonlinear system using FD Jacobian.

Solves the system:
    x1^2 + x2^2 = 1
    x1 - x2 = 0

Solutions: (1/√2, 1/√2) ≈ (0.7071, 0.7071) and (-1/√2, -1/√2)

Outputs:
- Table printed to stdout
- SVG plot of the system curves with solution marked
- Saved to test_output/demo_imsl_neqnf.svg
"""

from __future__ import annotations

import math
from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from nonlinearequations import solve_system_fd


def run_demo_imsl_neqnf() -> Dict[str, object]:
    """Run IMSL NEQNF example: 2-variable nonlinear system.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``solution`` (np.ndarray),
            ``fval`` (float), ``n_fev`` (int), and ``plot_path`` (str).
    """
    def f(x: np.ndarray) -> np.ndarray:
        return np.array([
            x[0]**2 + x[1]**2 - 1.0,  # unit circle
            x[0] - x[1],               # diagonal line
        ])

    x0 = np.array([0.8, 0.6])
    result = solve_system_fd(f, x0)
    sol = result.x
    residual = f(sol)
    expected = np.array([1 / math.sqrt(2), 1 / math.sqrt(2)])

    print("\nIMSL NEQNF Example: 2-Variable Nonlinear System")
    print("  System:  x1^2 + x2^2 = 1")
    print("           x1 - x2 = 0")
    print("-" * 55)
    print(f"{'Parameter':<25} {'Value':>20}")
    print("-" * 55)
    print(f"{'Starting point x0[0]':<25} {x0[0]:>20.6f}")
    print(f"{'Starting point x0[1]':<25} {x0[1]:>20.6f}")
    print(f"{'Solution x1':<25} {sol[0]:>20.10f}")
    print(f"{'Solution x2':<25} {sol[1]:>20.10f}")
    print(f"{'Expected (1/√2)':<25} {expected[0]:>20.10f}")
    print(f"{'Residual f1(x)':<25} {residual[0]:>20.2e}")
    print(f"{'Residual f2(x)':<25} {residual[1]:>20.2e}")
    print(f"{'|f(x)| (2-norm)':<25} {result.fval:>20.2e}")
    print(f"{'Function evals':<25} {result.n_fev:>20}")
    print(f"{'Converged':<25} {str(result.success):>20}")
    print("-" * 55)

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

    theta = np.linspace(0, 2 * np.pi, 400)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    diag_x = np.linspace(-1.2, 1.2, 200)
    diag_y = diag_x

    fig, ax = plt.subplots(figsize=(7, 6))
    ax.plot(circle_x, circle_y, color="#0e7490", linewidth=2.0, label=r"$x_1^2 + x_2^2 = 1$")
    ax.plot(diag_x, diag_y, color="#0891b2", linewidth=2.0, linestyle="--", label=r"$x_1 = x_2$")
    ax.scatter([sol[0], -sol[0]], [sol[1], -sol[1]], color="#d62728", s=100, zorder=5, label="Solutions")
    ax.annotate(
        f"({sol[0]:.4f}, {sol[1]:.4f})",
        xy=(sol[0], sol[1]),
        xytext=(sol[0] + 0.05, sol[1] - 0.15),
        fontsize=9, color="#d62728"
    )
    ax.annotate(
        f"({-sol[0]:.4f}, {-sol[1]:.4f})",
        xy=(-sol[0], -sol[1]),
        xytext=(-sol[0] - 0.35, -sol[1] + 0.1),
        fontsize=9, color="#d62728"
    )
    ax.set_aspect("equal")
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    ax.axhline(0, color="#aaaaaa", linewidth=0.7)
    ax.axvline(0, color="#aaaaaa", linewidth=0.7)
    ax.set_xlabel(r"$x_1$")
    ax.set_ylabel(r"$x_2$")
    ax.set_title("IMSL NEQNF: Nonlinear System Solve (FD Jacobian)")
    ax.legend()
    ax.grid(True, alpha=0.25)
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"solution": sol, "fval": result.fval, "n_fev": result.n_fev, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_neqnf()

Plot Output

The unit circle \(x_1^2 + x_2^2 = 1\) and the line \(x_1 = x_2\) intersect at the two solution points (red dots) identified by the solver.

Console Output

IMSL NEQNF Example: 2-Variable Nonlinear System
  System:  x1^2 + x2^2 = 1
           x1 - x2 = 0
-------------------------------------------------------
Parameter                                Value
-------------------------------------------------------
Starting point x0[0]                  0.800000
Starting point x0[1]                  0.600000
Solution x1                       0.7071067812
Solution x2                       0.7071067812
Expected (1/√2)                   0.7071067812
Residual f1(x)                        2.22e-16
Residual f2(x)                        0.00e+00
|f(x)| (2-norm)                       2.22e-16
Function evals                              10
Converged                                 True
-------------------------------------------------------