Solving a small nonlinear system with NumJac

A tiny worked example of using pymrm.NumJac to get a numerical Jacobian for a hand-written residual, and feeding the pair (g, jac) to a general-purpose solver. The system is

a\,x_0^2 + b\,x_1^2 - 1 = 0, \qquad x_0 - x_1 = 0,

(1)

solved once with scipy.optimize.root and once with pymrm.newton, to compare the two solvers on the same residual and Jacobian.

Imports¶

from pymrm import NumJac, newton
import scipy.optimize as optim
import numpy as np

Residual, Jacobian, and initial guess¶

numjac differentiates fun numerically; wrapper_fun adapts the (g, jac) pair to the dense Jacobian that scipy.optimize.root expects.

numjac = NumJac((2,))


def fun(x, a, b):
    return np.array([a * x[0] ** 2 + b * x[1] ** 2 - 1, x[0] - x[1]])


def wrapper_fun(x, *args):
    g, jac = numjac(lambda x: fun(x, *args), x)
    return g, jac.toarray()


x0 = np.array([0.5, 0.5])
args = (2, 3)

Solve with SciPy and with `pymrm.newton`¶

sol = optim.root(wrapper_fun, x0, args=args, jac=True)
print(f"Solution SciPy root:\n {sol}\n")

sol2 = newton(lambda x: numjac(lambda x: fun(x, *args), x), x0)
print(f"Solution pymrm newton:\n {sol2}\n")

Solution SciPy root:
  message: The solution converged.
 success: True
  status: 1
     fun: [ 2.220e-16  0.000e+00]
       x: [ 4.472e-01  4.472e-01]
  method: hybr
    nfev: 8
    njev: 1
    fjac: [[-8.769e-01 -4.807e-01]
           [ 4.807e-01 -8.769e-01]]
       r: [-2.080e+00 -1.841e+00  2.150e+00]
     qtf: [-3.263e-10  1.789e-10]

Solution pymrm newton:
    message: Converged
   success: True
       fun: [ 4.155e-10  0.000e+00]
         x: [ 4.472e-01  4.472e-01]
       nit: 4
 step_norm: 9.291678019311268e-11

Checking the result¶

With $x_0=x_1$ the first equation gives a closed form, $x_0=x_1=1/\sqrt{a+b}$ . Comparing both numerical solutions against it catches a wrong Jacobian, a solver that stopped short of convergence, or the wrong root of the (otherwise symmetric) system.

x_exact = 1.0 / np.sqrt(args[0] + args[1])
print(f"exact               x0=x1={x_exact:.12f}")
print(f"SciPy root          x0={sol.x[0]:.12f}, x1={sol.x[1]:.12f}, "
      f"max abs error {np.max(np.abs(sol.x - x_exact)):.3e}")
print(f"pymrm newton        x0={sol2.x[0]:.12f}, x1={sol2.x[1]:.12f}, "
      f"max abs error {np.max(np.abs(sol2.x - x_exact)):.3e}")

exact               x0=x1=0.447213595500
SciPy root          x0=0.447213595500, x1=0.447213595500, max abs error 5.551e-17
pymrm newton        x0=0.447213595500, x1=0.447213595500, max abs error 1.110e-16

Imports¶

Residual, Jacobian, and initial guess¶

Solve with SciPy and with pymrm.newton¶

Checking the result¶

Solve with SciPy and with `pymrm.newton`¶