Solve an equation using a python numerical solver in numpy

Question

I have an equation, as follows:

R - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau))) = 0.

I want to solve for tau in this equation using a numerical solver available within numpy. What is the best way to go about this?

The values for R and a in this equation vary for different implementations of this formula, but are fixed at particular values when it is to be solved for tau.

Answer 1

In conventional mathematical notation, your equation is

The SciPy fsolve function searches for a point at which a given expression equals zero (a "zero" or "root" of the expression). You'll need to provide fsolve with an initial guess that's "near" your desired solution. A good way to find such an initial guess is to just plot the expression and look for the zero crossing.

#!/usr/bin/python

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

# Define the expression whose roots we want to find

a = 0.5
R = 1.6

func = lambda tau : R - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau))) 

# Plot it

tau = np.linspace(-0.5, 1.5, 201)

plt.plot(tau, func(tau))
plt.xlabel("tau")
plt.ylabel("expression value")
plt.grid()
plt.show()

# Use the numerical solver to find the roots

tau_initial_guess = 0.5
tau_solution = fsolve(func, tau_initial_guess)

print "The solution is tau = %f" % tau_solution
print "at which the value of the expression is %f" % func(tau_solution)

Answer 2

You can rewrite the equation as

• For integer a and non-zero R you will get a solutions in the complex space;
• There are analytical solutions for a=0,1,...4(see here);

So in general you may have one, multiple or no solution and some or all of them may be complex values. You may easily throw scipy.root at this equation, but no numerical method will guarantee to find all the solutions.

To solve in the complex space:

import numpy as np
from scipy.optimize import root

def poly(xs, R, a):
    x = complex(*xs)
    err = R * x - x + 1 - R
    return [err.real, err.imag]

root(poly, x0=[0, 0], args=(1.2, 6))