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.
In Python, NumPy (Numerical Python), SciPy (Scientific Python) and SymPy (Symbolic Python) libraries can be used to solve systems of linear equations.
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)
You can rewrite the equation as
a
and non-zero R
you will get a
solutions in the complex space;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))
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With