Search code examples
pythonscipyequationsolverroot-finding

How to solve a polynomial equation for one unknown in python?


I have been trying to solve the following equation for R1 in python using scipy, but I get the error:

result = root_scalar(equation, bracket=[0, 8])
Traceback (most recent call last):

  Cell In[1638], line 1
    result = root_scalar(equation, bracket=[0, 8])

  File ~\anaconda3\lib\site-packages\scipy\optimize\_root_scalar.py:275 in root_scalar
    r, sol = methodc(f, a, b, args=args, **kwargs)

  File ~\anaconda3\lib\site-packages\scipy\optimize\_zeros_py.py:784 in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)

ValueError: f(a) and f(b) must have different signs

The equation is:

K= 10**(a0 + (a1 * np.log10(R1/R2)) + (a2 * (np.log10(R1/R2))**2) + (a3 * (np.log10(R1/R2))**3) + (a4 * (np.log10(R1/R2))**4)) + 0.0166

The code that I have been using is the following:

R2 = 0.002294000005349517
K = 0.09539999812841415

a0 = -1.1358
a1 = -2.1146
a2 = 1.6474
a3 = -1.1428
a4 = -0.6190

def equation(log_R1_R2):
    x = log_R1_R2
    return 10**(a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4) - (K - 0.0166)

result = root_scalar(equation, bracket=[0, 8])
log_R1_R2 = result.root
R1 = R2 * 10**log_R1_R2

Solution

  • This problem seems better suited for a polynomial root solver. Using the numpy.polynomial.Polynomial class, we can create the polynomial in question, get the roots, and find R1 using those roots. As for which roots you want to use, that's for you to decide since I don't know the application.

    import numpy as np
    from numpy.polynomial import Polynomial
    
    R2 = 0.002294000005349517
    K = 0.09539999812841415
    
    a0 = -1.1358
    a1 = -2.1146
    a2 = 1.6474
    a3 = -1.1428
    a4 = -0.6190
    
    p = Polynomial([-np.log10(K - 1/60) + a0, a1, a2, a3, a4])
    roots = p.roots()
    print(roots)
    R1 = R2*10**roots
    

    Output (i.e. the roots):

    [-3.07249423+0.j, -0.0149377 +0.j, 0.62061419-0.86009349j, 0.62061419+0.86009349j]