9c. Numerical one-dimensional bounded root finding

---

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
# ------------------------------

import numpy as np
import scipy.optimize

---

When finding the fixed point of the repressilator system that includes mRNA in Chapter 9, we were tasked with finding the value of \(x_0\) that is the solution to

\begin{align} (x_0 - \beta\rho)(1+x_0^n) = \beta. \end{align}

Finding an analytical solution of this equation for \(x_0\) is usually impossible. (It is possible when \(n\) is an integer between zero and three.) We therefore seek a numerical method for finding \(x_0\). This is an example of a root finding problem. It can be cast into a problem of finding \(f(x) = 0\) for some function \(f(x)\). Do not confuse this \(f(x)\) with regulatory functions we have used; we are using \(f(x)\) here to be an arbitrary function. In the present case, \(f(x) = \beta - (x_0 - \beta\rho)(1+x_0^n)\).

We will explore algorithms for more general multidimensional root finding problem in later chapters, but for now, we seek a scalar \(x_0\). In this case, we know a lot about the fixed point. We know from our analysis in Chapter 9 that it exists and is unique. We also know that it lies between \(x_0 = 0\) (where \(f(0) = \beta(1+\rho) > 0\)) and \(x_0 = \beta(1+\rho)\) (since \(f(\beta(1+\rho)) < 0\)). When we have bounds and guarantees of uniqueness for a scalar root, we can use a bisection method to find the root. The bisection method is guaranteed to find the root of the function \(f(x)\) on an interval \([a,b]\), provided \(f(x)\) is continuous and \(f(a)\) and \(f(b)\) have opposite sign, which is the case here.

Brent’s method also has this guarantee, but is more efficient that using bisection. We will not go into the details of these respective algorithms, and leave you to read the Wikipedia pages we have linked to, which provide quite clear descriptions.

Brent’s method is available in the scipy.optimize.brentq() function. It takes as an argument the function \(f(x)\) whose root is to be found, and the left and right bounds for the root. We can write a function to find the fixed point for given \(\beta\), \(\rho\), and \(n\).

def fixed_point(beta, n, rho):
    def _root_function(x):
        return beta - (x - beta * rho) * (1 + x ** n)

    return scipy.optimize.brentq(_root_function, 0, beta * (1 + rho))

As an example, we can put it to use for values of \(\beta\), \(n\), and \(\rho\) for which an analytical solution is not possible.

beta = 10.07
n = 4.78
rho = 0.34

fixed_point(beta, n, rho)

3.4507547941440446

Computing environment

%load_ext watermark
%watermark -v -p numpy,scipy,jupyterlab

Python implementation: CPython
Python version       : 3.10.10
IPython version      : 8.12.0

numpy     : 1.23.5
scipy     : 1.10.1
jupyterlab: 3.5.3