Multiphase Reactor Modeling using PyMRM

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


# Define the objective function and its Jacobian
numjac = NumJac((2,))

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

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

# Solve the system
def wrapper_fun(x, *args):
    g, Jac = numjac(lambda x: fun(x, *args), x)
    return g, Jac.toarray()

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