Python:Finding roots

Introduction

The SciPy package scipy.optimize has several routines for finding roots of equations. This document is going to focus on the brentq function for finding the root of a single-variable continuous function. The function can only find one root at a time and it requires brackets for the root. For open root-finding, use fsolve. All the topions below for brentq work with fsolve, the only difference is instead of giving a left and right bracket you give a single initial guess.

All of the examples assume that

import numpy as np
import scipy.optimization as opt

have been called.

Most General Case

The following examples show a method that will work regardless of how many input variables your function has or how you define your function - whether it is a built-in or user-defined function. The paradigm is:

ROOT = opt.brentq(lambda DUMMY_VAR: FUNCTION_THING, LEFT_BRACKET, RIGHT_BRACKET)

where

• ROOT is the calculated value of the requested variable when the function is 0
• DUMMY_VAR is the variable you want to use in this FUNCTION_THING to indicate which of the various inputs brentq is allowed to alter
• LEFT_BRACKET and RIGHT_BRACKET are values of DUMMY that bracket the root - that is, where FUNCTION_THING has different signs
  • Note that brentq can only find one root at a time, so you cannot load this brackets with several values and try to have brentq find multiple results; for that, you need loops
• FUNCTION_THING can be a built-in function, a user-defined function, or a calculation - whatever it is that you are trying to set equal to zero; note that DUMMY_VAR must appear somewhere in this expression for brentq to be able to do anything.

Example

If you want to solve \(f(x)=(x+4)(x-1)=x^2+3x-4=0\) for instance, you have several options:

Function on the Fly

If you only want to solve this problem a couple times, and the calculation only requires one line of code, the process with the least "overhead" involves creating the function on the fly inside the brentq command. All you have to do is put your expression in for the FUNCTION_THING, making sure the DUMMY_VAR is appropriately used - you can use a lambda function for that:

r1 = opt.brentq(lambda x: x**2+3*x-4, -5, -3)
r2 = opt.brentq(lambda x: x**2+3*x-4, 0, 3)

where -5 and -3 were chosen since they bracket the root at -4 and 0 and 3 were chosen since they bracket the root at 1.

Defined Functions

If you are going to solve roots multiple times, you may want to define a function and then call it using a lambda function:

def fun(x):
    return x**2+3*x-4

r1 = opt.brentq(lambda xi: fun(xi), -4, -2)
r2 = opt.brentq(lambda xi: fun(xi), 0, 3)

Note that the dummy variable xi could have been an x here.

Single-Variable Functions

If you have a function that takes a single variable, there is a slightly simpler way to run brentq - all you need to do is give the name of the variable to FUNCTION_THING rather than setting up the DUMMY_VAR:

def fun(x):
    return x**2+3*x-4

r1 = opt.brentq(fun, -4, -2)
r2 = opt.brentq(fun, 0, 3)

Multiple-Variable Functions

If you have a function with multiple variables, first remember that all the other variables need to have known values. For example, if you have \(f(x,y)=\cos(x)+\sin(y)\), and you want to know where \(f(x,y)=0\), you need to know either \(x\) or \(y\).

Unknown as first parameter in function call

Imagine you want to know \(x\) when \(y\) is 2.5 radians - you can solve this the following ways:

On the fly

xval = opt.brentq(lambda x: np.cos(x) + np.sin(2.5), 0, np.pi)

Defined function

def fun2(x,y):
    return np.cos(x)+np.sin(y)

xval = opt.brentq(lambda x: fun2(x, 2.5), 0, np.pi)

or, since the variable you are solving for is first, you can give brentq a tuple with the 2nd (and beyond) variable values:

def fun2(x,y):
    return np.cos(x)+np.sin(y)

xval = opt.brentq(fun2, 0, np.pi, args=(2.5,)) # formally
xval = opt.brentq(fun2, 0, np.pi, 2.5) # less formally since the args positionally come after the right bracket

Unknown not first parameter in function call

Now imagine you want to know \(y\) when \(x\) is 2 - you have to create a lambda function for that, either on the fly or based on the defined function:

On the fly

yval = opt.brentq(lambda x: np.cos(2) + np.sin(y), 0, 1)

= Defined Function

def fun2(x,y):
    return np.cos(x)+np.sin(y)

yval = opt.brentq(lambda y: fun2(2, y), 0, 1)

Questions

Post your questions by editing the discussion page of this article. Edit the page, then scroll to the bottom and add a question by putting in the characters *{{Q}}, followed by your question and finally your signature (with four tildes, i.e. ~~~~). Using the {{Q}} will automatically put the page in the category of pages with questions - other editors hoping to help out can then go to that category page to see where the questions are. See the page for Template:Q for details and examples.

References