NEQNJ Example — Nonlinear System with Analytic Jacobian¶
This example solves a 2-variable nonlinear system using an analytic Jacobian (IMSL NEQNJ):
\[ \begin{align}\begin{aligned}f_1(x, y) = x^2 + y^2 - 4 = 0 \quad \text{(circle of radius 2)}\\f_2(x, y) = xy - 1 = 0 \quad \text{(rectangular hyperbola)}\end{aligned}\end{align} \]
The analytic Jacobian:
\[\begin{split}J = \begin{pmatrix} 2x & 2y \\ y & x \end{pmatrix}\end{split}\]
is supplied directly to nonlinearequations.solve_system(), reducing
function evaluations compared to the finite-difference variant. Starting
from \((1.5, 0.5)\) and \((-1.5, -0.5)\), the solver converges to
both intersection points near \((\approx 1.9319,\, 0.5176)\) and its
symmetric counterpart.
Example Code¶
"""IMSL NEQNJ example: solve a nonlinear system with an analytic Jacobian.
Solves the system:
f1(x, y) = x^2 + y^2 - 4 (circle of radius 2)
f2(x, y) = x*y - 1
Two solutions near (1.93, 0.52) and (-1.93, -0.52).
Uses :func:`nonlinearequations.solve_system` (IMSL NEQNJ) which accepts an
explicit analytic Jacobian, reducing function evaluations versus the
finite-difference variant.
Outputs:
- Table printed to stdout
- SVG contour plot of f1=0 and f2=0 with solution marked
- Saved to test_output/example_imsl_neqnj.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 solve_system
def run_demo_imsl_neqnj() -> Dict[str, object]:
"""Run IMSL NEQNJ example: 2-variable nonlinear system with analytic Jacobian.
Solves f1 = x^2 + y^2 - 4 = 0 and f2 = x*y - 1 = 0 from two different
starting points to recover both solution branches.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``solutions`` (list of
np.ndarray), ``fvals`` (list of float), ``n_fev`` (int), and
``plot_path`` (str).
"""
def f(v: np.ndarray) -> np.ndarray:
"""System function for the nonlinear equations.
Args:
v (np.ndarray): Current point [x, y].
Returns:
np.ndarray: Residuals [f1, f2].
"""
x, y = v[0], v[1]
return np.array([
x**2 + y**2 - 4.0, # circle of radius 2
x * y - 1.0, # rectangular hyperbola
])
def jac(v: np.ndarray) -> np.ndarray:
"""Analytic Jacobian of the system.
Args:
v (np.ndarray): Current point [x, y].
Returns:
np.ndarray: 2x2 Jacobian matrix df_i/dx_j.
"""
x, y = v[0], v[1]
return np.array([
[2.0 * x, 2.0 * y], # d(f1)/dx, d(f1)/dy
[y, x], # d(f2)/dx, d(f2)/dy
])
# Solve from two starting points to get both solutions
x0_pos = np.array([1.5, 0.5])
x0_neg = np.array([-1.5, -0.5])
result_pos = solve_system(f, jac, x0_pos)
result_neg = solve_system(f, jac, x0_neg)
solutions = [result_pos.x, result_neg.x]
fvals = [result_pos.fval, result_neg.fval]
total_fev = result_pos.n_fev + result_neg.n_fev
print("\nIMSL NEQNJ Example: Nonlinear System with Analytic Jacobian")
print(" System: f1: x^2 + y^2 = 4 (circle, radius 2)")
print(" f2: x*y = 1 (hyperbola)")
print("-" * 60)
print(f"{'Parameter':<30} {'Value':>20}")
print("-" * 60)
for i, (sol, fv) in enumerate(zip(solutions, fvals), start=1):
print(f"{'Solution ' + str(i) + ' x':<30} {sol[0]:>20.10f}")
print(f"{'Solution ' + str(i) + ' y':<30} {sol[1]:>20.10f}")
print(f"{'Residual norm ' + str(i):<30} {fv:>20.2e}")
print(f"{'Total function evals':<30} {total_fev:>20}")
print(f"{'Both converged':<30} {str(result_pos.success and result_neg.success):>20}")
print("-" * 60)
# ------------------------------------------------------------------
# Plot: contour lines f1=0 and f2=0 with solutions marked
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_neqnj.svg"
x_grid = np.linspace(-2.8, 2.8, 500)
y_grid = np.linspace(-2.8, 2.8, 500)
X, Y = np.meshgrid(x_grid, y_grid)
F1 = X**2 + Y**2 - 4.0
F2 = X * Y - 1.0
fig, ax = plt.subplots(figsize=(7, 6))
ax.contour(X, Y, F1, levels=[0], colors=["#0e7490"], linewidths=2.0)
ax.contour(X, Y, F2, levels=[0], colors=["#0891b2"], linewidths=2.0, linestyles="--")
# Dummy lines for legend
ax.plot([], [], color="#0e7490", linewidth=2.0, label=r"$f_1 = x^2 + y^2 - 4 = 0$")
ax.plot([], [], color="#0891b2", linewidth=2.0, linestyle="--",
label=r"$f_2 = xy - 1 = 0$")
sol_x = [s[0] for s in solutions]
sol_y = [s[1] for s in solutions]
ax.scatter(sol_x, sol_y, color="#d62728", s=100, zorder=5, label="Solutions (NEQNJ)")
for sol in solutions:
ax.annotate(
f"({sol[0]:.4f}, {sol[1]:.4f})",
xy=(sol[0], sol[1]),
xytext=(sol[0] + 0.12, sol[1] - 0.25),
fontsize=9,
color="#d62728",
)
ax.axhline(0, color="#aaaaaa", linewidth=0.7)
ax.axvline(0, color="#aaaaaa", linewidth=0.7)
ax.set_xlim(-2.8, 2.8)
ax.set_ylim(-2.8, 2.8)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("IMSL NEQNJ: Nonlinear System with Analytic Jacobian")
ax.legend(loc="upper left")
ax.grid(True, alpha=0.25)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"solutions": solutions,
"fvals": fvals,
"n_fev": total_fev,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_neqnj()
Plot Output¶
Contour lines where \(f_1 = 0\) (circle, teal) and \(f_2 = 0\) (hyperbola, dashed) intersect at the two solution points (red dots) found by the analytic-Jacobian solver.¶
Console Output¶
IMSL NEQNJ Example: Nonlinear System with Analytic Jacobian
System: f1: x^2 + y^2 = 4 (circle, radius 2)
f2: x*y = 1 (hyperbola)
------------------------------------------------------------
Parameter Value
------------------------------------------------------------
Solution 1 x 1.9318516526
Solution 1 y 0.5176380902
Residual norm 1 6.99e-13
Solution 2 x -1.9318516526
Solution 2 y -0.5176380902
Residual norm 2 6.99e-13
Total function evals 22
Both converged True
------------------------------------------------------------