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.
    """
    # 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 "Invalid input: equations must be a non-empty 2D list of strings."

    # Flatten equations
    flat_equations = []
    for row in equations:
        if not isinstance(row, list):
            return "Invalid input: 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 "Invalid input: equations must contain at least one string."

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

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

    try:
        x0 = np.array([float(v) for v in initial_guess], dtype=float)
    except (TypeError, ValueError):
        return "Invalid input: x_zero must contain numeric values."

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

    if n_vars != n_eqs:
        return f"Invalid input: 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
    try:
        x = ca.SX.sym('x', n_vars)
    except Exception as exc:
        return f"Error: failed to create symbolic variable: {exc}"

    # Parse equations as CasADi expressions
    try:
        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
            }
            try:
                # Evaluate expression with CasADi variables
                expr_result = eval(eq_str, {"__builtins__": {}, **eval_context}, {})
                residuals.append(expr_result)
            except Exception as exc:
                return f"Error: failed to parse equation {i}: {exc}"
    except Exception as exc:
        return f"Error: equation parsing failed: {exc}"

    # Set up optimization problem to minimize sum of squares of residuals
    try:
        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:
                        try:
                            lbx[i] = float(bound[0])
                            ubx[i] = float(bound[1])
                        except (TypeError, ValueError):
                            return f"Invalid input: bounds for variable {i} must be numeric or None."

        # 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 exc:
        return f"Error: {exc}"

Overview

Excel Usage

Examples

Python Code

Online Calculator