Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Solve an equation using a python numerical solver in numpy

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.

like image 626
stars83clouds Avatar asked Mar 30 '14 10:03

stars83clouds


People also ask

Can Python solve systems of equations?

In Python, NumPy (Numerical Python), SciPy (Scientific Python) and SymPy (Symbolic Python) libraries can be used to solve systems of linear equations.


2 Answers

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) 
like image 56
nibot Avatar answered Sep 22 '22 10:09

nibot


You can rewrite the equation as

eq

  • 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)) 
like image 35
behzad.nouri Avatar answered Sep 23 '22 10:09

behzad.nouri