Solve an equation using a python numerical solver in numpy

stars83clouds picture stars83clouds · Mar 30, 2014 · Viewed 84.9k times · Source

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

nibot picture nibot · Mar 30, 2014

In conventional mathematical notation, your equation is

$$ R = \frac{1 - e^{-\tau}}{1 - e^{-a\cdot\tau}}$$

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)