How do I find the roots of the function $f(\beta) = \gamma + [1-e^{-j\beta}]/[1-e^{(-j+1)\beta}]$, using python

Question

I have the following equation which I cannot find a closed form solution for $\beta$, thus I would like to solve for $\beta$ computationally:

$$\gamma = \frac{1-e^{-jT\beta}}{1-e^{-(j+1)T\beta}}$$

variables $\gamma$, $j$ and $T$ are known. Unfortunately I do not have knowledge on the brackets of the roots of the equation.

Is there an algorithm in python which is able to solve for the roots without knowledge on the brackets?

Answer 1

Scipy has several options for root finding algorithms: https://docs.scipy.org/doc/scipy-0.13.0/reference/optimize.html#root-finding

These algorithms tend to want either bounds for the possible solution or the derivative. In this case, bounds seemed like less work :)

edit: After re-reading your question, It sounds like you don't have an estimate for the upper and lower bounds. In that case, I suggest Newton's method (scipy.optim.newton, in which case you'll have to compute the derivative (by hand or with sympy).

import matplotlib
matplotlib.use('Agg')

import numpy as np
from scipy.optimize import brentq
import matplotlib.pyplot as plt

def gamma_func(beta, j, T):
    return (1-np.exp(-j*T*beta)) / (1-np.exp(-(j+1)*T*beta))

# Subtract gamma from both sides so that we have f(beta,...) = 0.
# Also, set beta, the independent variable, to the first
# argument for compatibility with solver
def f(beta, gamma, j, T):
    return gamma_func(beta, j, T) - gamma


# Parameters
gamma = 0.5
j = 2.3
T = 1.5

# Lower and upper bounds for solution
beta_low = -3
beta_high = 3

# Error tolerance for solution
tol = 1e-4

# True solution
beta_star = brentq(f, beta_low, beta_high, args=(gamma, j, T), xtol=tol)


# Plot the solution
plt.figure()
plt.clf()
beta_arr = np.linspace(-10, 10, 1001)
gamma_arr = gamma_func(beta_arr, j, T)
plt.plot(beta_arr, gamma_arr)
plt.plot(beta_star, gamma_func(beta_star, j, T), 'o')
plt.savefig('gamma_plot.png')

# Print the solution
print("beta_star = {:.3f}".format(beta_star))
print("gamma(beta_star) = {:.3f}".format(gamma_func(beta_star, j, T)))

beta_star = -0.357
gamma(beta_star) = 0.500

See the documentation for the other solvers, including multidimensional solvers.

Good luck!

Oliver