SciPy - Finding Roots of Multivariate Functions

Finding Roots of Multivariate Functions in SciPy refers to the process of finding the values of multiple variables, typically represented as a vector that satisfy a system of nonlinear equations i.e., solving for the root or zero of a multivariate function. In other words we can say it involves finding the points where a vector-valued function F(x)=0 where F(x) is a set of equations and x a vector of variables.

How SciPy Helps in Finding Roots?

In SciPy we have the function scipy.optimize.root() which is used to find the roots of systems of multivariate nonlinear equations. It is a part of the scipy.optimize module which provides optimization and root-finding routines. The root function allows us to solve these systems by providing an initial guess and choosing from a variety of methods for solving the system.

Why to Use SciPy to find Multivariate Root?

Here are the reasons why to use SciPy for Multivariate Root Finding −

• Robust and Flexible Algorithms: SciPy provides various solvers such as 'hybr', 'lm', 'broyden1', etc which tailored for different nonlinear systems by making it versatile for diverse problems.
• Handles Complex Systems: It can solve coupled, nonlinear equations that lack analytical solutions.
• Ease of Use: A simple interface lets users define equations, provide initial guesses, and solve efficiently.
• Automatic or Custom Jacobians: SciPy can approximate derivatives or accept user-provided Jacobians for faster convergence.
• Scalable for Large Systems: Solvers such as 'broyden1' handle large systems without needing explicit Jacobians.
• Reliable and Optimized: Built on trusted libraries like MINPACK by ensuring accuracy and performance.
• Integration: Works seamlessly with NumPy, Matplotlib and SymPy for end-to-end workflows.
• Diagnostic Tools: Provides detailed outputs for convergence, residuals, and solution validation.

Syntax

Following is the Syntax of the function scipy.optimize.root() which is used to find the roots of multivariate functions −

scipy.optimize.root(
   fun, 
   x0, 
   args=(), 
   method='hybr', 
   jac=None, 
   tol=None, 
   options={}
)

Parameters

• fun(callable): The function whose roots are sought. It should take an array x and return an array of the same shape.
• x0(array_like): Initial guess for the solution.
• args(tuple, optional): Extra arguments to pass to the function fun.
• method(str, optional): The algorithm to use. Options include 'hybr', 'lm', 'broyden1', 'anderson', etc.
• jac(callable, optional): A function to compute the Jacobian matrix of fun. If not provided then finite differences are used.
• tol(float, optional): Tolerance for termination.
• options(dict, optional): Solver-specific options like max iterations or verbosity level.

Return Value

scipy.optimize.root() returns an OptimizeResult object containing the below values −

• x: The solution array.
• success: A boolean indicating if the solver succeeded.
• message: Description of the cause of termination.
• fun: The residuals at the solution.

Finding the Root Multiple functions

Here in this example we will solve for the roots of multiple function involving multiple variables for the equations f₁(x,y) = e^x+y²-1 = 0 and f₂(x,y) = x³+3x²+1 = 0 −

import numpy as np
from scipy.optimize import root

# Define the system of nonlinear equations
def fun_system(vars):
    x, y = vars
    # f1(x, y) = e^x + y^2 - 1
    # f2(x, y) = x^3 - 3x + 1
    f1 = np.exp(x) + y**2 - 1
    f2 = x**3 - 3*x + 1
    return [f1, f2]  # Return both residuals

# Initial guess for [x, y]
initial_guess = [0.0, 0.0]

# Use root to solve the system of equations
result = root(fun_system, initial_guess)

# Display the result
if result.success:
    print(f"Roots found: x = {result.x[0]}, y = {result.x[1]}")
else:
    print(f"Root finding failed: {result.message}")

Here is the output of the function scipy.optimize.root() which is used to find the root of multiple functions −

Root finding failed: The iteration is not making good progress, as measured by the
  improvement from the last ten iterations.

Solving a System with More Equations & Variables

Here is an example which shows how to solve a system with more equations and variables with the help of the function scipy.optimize.root() −

import numpy as np
from scipy.optimize import root

# Define the system of equations
def system_of_equations(x):
   return [
      x[0]**2 + x[1]**2 + x[2]**2 - 9,  # x + x + x = 9
      x[0] + x[1] - x[2] - 1,            # x + x - x = 1
      x[0] - x[1] - x[2] - 2             # x - x - x = 2
   ]

# Initial guess for the solution
x0 = [1, 1, 1]

# Solve the system using the 'hybr' method (default)
result = root(system_of_equations, x0, method='hybr')

# Check the result and display the solution
if result.success:
    print("Solution found:", result.x)
else:
    print("Solver failed:", result.message)

Following is the output of the function scipy.optimize.root() which is used to solve a system with more equations and variables −

Solution found: [ 2.70256242 -0.5         1.20256242]