CA_ROOT

Overview

The CA_ROOT function solves systems of nonlinear equations using CasADi, an open-source symbolic framework for numerical optimization. Given a set of equations f_1(x) = 0, f_2(x) = 0, \ldots, f_n(x) = 0, the function finds values of x = [x_0, x_1, \ldots, x_{n-1}] that simultaneously satisfy all equations.

Root finding is the process of determining where a function equals zero. For systems of nonlinear equations, this is typically solved by reformulating the problem as minimizing the sum of squared residuals:

\min_{x} \sum_{i=1}^{n} f_i(x)^2

This implementation uses CasADi’s nlpsol interface with the Sequential Quadratic Programming (SQP) method (sqpmethod) as the default solver. CasADi automatically computes the Jacobian matrix using algorithmic differentiation, which provides exact derivative information without the numerical errors associated with finite differences. For details, see the CasADi documentation and the CasADi GitHub repository.

The function supports:

Arbitrary system size: Solve systems with any number of equations and variables (must be equal)
Optional bounds: Constrain solutions to specific regions using lower and upper bounds on each variable
Custom solvers: Choose alternative NLP solvers via the solver parameter

Equations are specified as strings using the notation x[0], x[1], etc. for variables, with support for standard mathematical operators and functions including sqrt, sin, cos, exp, log, and exponentiation via ^ or **.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=CA_ROOT(equations, x_zero, bounds, solver)

equations (list[list], required): 2D list of equation strings in terms of x[0], x[1], etc. Each expression should equal zero at the solution.
x_zero (list[list], required): 2D list containing initial guess values for each variable. Must have same length as number of equations.
bounds (list[list], optional, default: null): 2D list of [lower, upper] bound pairs for each variable to constrain the solution space.
solver (str, optional, default: “sqpmethod”): Solver name to use for optimization.

Returns (list[list]): 2D list [[x1, x2, …]], or error message string.

Examples

Example 1: Circle and line intersection

Inputs:

equations	x_zero
x[0]^2 + x[1]^2 - 1	1	1
x[0] - x[1]

Excel formula:

=CA_ROOT({"x[0]^2 + x[1]^2 - 1";"x[0] - x[1]"}, {1,1})

Expected output:

Result
0.7071	0.7071

Example 2: Coupled cubic system

Inputs:

equations	x_zero
x[0]^3 + x[1] - 1	0.7	0.7
x[0] + x[1]^3 - 1

Excel formula:

=CA_ROOT({"x[0]^3 + x[1] - 1";"x[0] + x[1]^3 - 1"}, {0.7,0.7})

Expected output:

Result
0.6826	0.6826

Example 3: Circle and line with bounds

Inputs:

equations	x_zero		bounds
x[0]^2 + x[1]^2 - 1	-0.5	-0.5	-1	0
x[0] - x[1]			-1	0

Excel formula:

=CA_ROOT({"x[0]^2 + x[1]^2 - 1";"x[0] - x[1]"}, {-0.5,-0.5}, {-1,0;-1,0})

Expected output:

Result
-0.7071	-0.7071

Example 4: Three-variable system with all arguments

Inputs:

equations	x_zero				solver
x[0]^2 + x[1]^2 + x[2]^2 - 1	0.5	0.5	0.5	1	sqpmethod
x[0] - x[1]				1
x[1] - x[2]				1

Excel formula:

=CA_ROOT({"x[0]^2 + x[1]^2 + x[2]^2 - 1";"x[0] - x[1]";"x[1] - x[2]"}, {0.5,0.5,0.5}, {0,1;0,1;0,1}, "sqpmethod")

Expected output:

Result
0.5773	0.5773	0.5773

Python Code

import re
import math
import casadi as ca
import numpy as np

def ca_root(equations, x_zero, bounds=None, solver='sqpmethod'):
    """
    Solve a system of nonlinear equations using CasADi with automatic Jacobian.

    See: https://web.casadi.org/docs/

    This example function is provided as-is without any representation of accuracy.

    Args:
        equations (list[list]): 2D list of equation strings in terms of x[0], x[1], etc. Each expression should equal zero at the solution.
        x_zero (list[list]): 2D list containing initial guess values for each variable. Must have same length as number of equations.
        bounds (list[list], optional): 2D list of [lower, upper] bound pairs for each variable to constrain the solution space. Default is None.
        solver (str, optional): Solver name to use for optimization. Default is 'sqpmethod'.

    Returns:
        list[list]: 2D list [[x1, x2, ...]], or error message string.
    """
    try:
      # Normalize inputs to 2D list
      def to2d(x):
          return [[x]] if not isinstance(x, list) else x

      equations = to2d(equations)
      x_zero = to2d(x_zero)

      # Validate equations input
      if not isinstance(equations, list) or len(equations) == 0:
          return "Error: equations must be a non-empty 2D list of strings."

      # Flatten equations
      flat_equations = []
      for row in equations:
          if not isinstance(row, list):
              return "Error: equations must be a 2D list of strings."
          for eq in row:
              if isinstance(eq, str) and eq.strip():
                  flat_equations.append(eq)

      if len(flat_equations) == 0:
          return "Error: equations must contain at least one string."

      # Validate and extract initial guess
      if not isinstance(x_zero, list) or len(x_zero) == 0:
          return "Error: x_zero must be a 2D list."

      initial_guess = x_zero[0]
      if not isinstance(initial_guess, list) or len(initial_guess) == 0:
          return "Error: x_zero must contain at least one value."

      x0 = np.array([float(v) for v in initial_guess], dtype=float)

      n_vars = len(x0)
      n_eqs = len(flat_equations)

      if n_vars != n_eqs:
          return f"Error: number of variables ({n_vars}) must match equations ({n_eqs})."

      # Convert caret to double asterisk for exponentiation
      equations_py = [re.sub(r'\^', '**', eq) for eq in flat_equations]

      # Create CasADi symbolic variable
      x = ca.SX.sym('x', n_vars)

      # Parse equations as CasADi expressions
      residuals = []
      for i, eq_str in enumerate(equations_py):
          # Build evaluation context with CasADi functions
          eval_context = {
              "x": x,
              "sqrt": ca.sqrt,
              "sin": ca.sin,
              "cos": ca.cos,
              "tan": ca.tan,
              "exp": ca.exp,
              "log": ca.log,
              "log10": ca.log10,
              "abs": ca.fabs,
              "fabs": ca.fabs,
              "asin": ca.asin,
              "acos": ca.acos,
              "atan": ca.atan,
              "sinh": ca.sinh,
              "cosh": ca.cosh,
              "tanh": ca.tanh,
              "pow": ca.power,
              "pi": ca.pi,
              "e": math.e
          }
          # Evaluate expression with CasADi variables
          expr_result = eval(eq_str, {"__builtins__": {}, **eval_context}, {})
          residuals.append(expr_result)

      F = ca.Function('F', [x], residuals)
      residual_vector = ca.vertcat(*residuals)
      objective = ca.dot(residual_vector, residual_vector)

      # Set up NLP
      nlp = {'x': x, 'f': objective}

      # Add bounds if provided
      lbx = -np.inf * np.ones(n_vars)
      ubx = np.inf * np.ones(n_vars)

      if bounds is not None:
          if isinstance(bounds, list) and len(bounds) == n_vars:
              for i, bound in enumerate(bounds):
                  if isinstance(bound, list) and len(bound) >= 2:
                      lbx[i] = float(bound[0])
                      ubx[i] = float(bound[1])

      # Create and configure solver
      opts = {
          'qpsol': 'qrqp',
          'qpsol_options': {'print_iter': False, 'print_header': False},
          'print_time': 0,
          'max_iter': 3000
      }

      solver_obj = ca.nlpsol('solver', solver, nlp, opts)

      # Solve
      result = solver_obj(x0=x0, lbx=lbx, ubx=ubx)

      if result is None:
          return "Error: solver returned None."

      x_opt = result['x'].full().flatten()

      # Verify solution quality
      residuals_at_opt = F(x_opt)
      # Handle single vs multiple residuals - CasADi returns DM for single output
      if n_eqs == 1:
          residual_norm = float(ca.fabs(residuals_at_opt))
      else:
          residual_norm = float(ca.norm_2(ca.vertcat(*residuals_at_opt)))

      if residual_norm > 1e-3:
          return f"Error: solver did not converge (residual norm: {residual_norm:.6f})."

      solution = [float(val) for val in x_opt]
      return [solution]
    except Exception as e:
      return f"Error: {str(e)}"

Overview

Excel Usage

Examples

Python Code

Online Calculator