# Using numpy polynomial module - is there a better way?

I'm working on a project where I need to solve for one of the roots of a quartic polymonial many, many times. Is there a better, i.e., faster way to do this? Should I write my own C-library? The example code is below.

# this code calculates the pH of a solution as it is
# titrated with base and then plots it.

import numpy.polynomial.polynomial as poly
import numpy as np
import matplotlib.pyplot as plt

# my pH calculation function
# assume two distinct pKa's  solution is a quartic equation
def pH(base, Facid1, Facid2):
ka1 = 2.479496e-6
ka2 = 1.87438e-9
kw = 1.019230e-14
a = 1
b = ka1+ka2+base
c = base*(ka1+ka2)-(ka1*Facid1+ka2*Facid2)+ka1*ka2-kw
d = ka1*ka2*(base-Facid1-Facid2)-kw*(ka1+ka2)
e = -kw*ka1*ka2
p = poly.Polynomial((e,d,c,b,a))
return -np.log10(p.roots()[3]) #only need the 4th root here

# Define the concentration parameters
Facid1 = 0.002
Facid2 = 0.001
Fbase = 0.005    #the maximum base addition

# Generate my vectors
x = np.linspace(0., Fbase, 200)
y = [pH(base, Facid1, Facid2) for base in x]

# Make the plot frame
fig = plt.figure()

# Set the limits
ax.set_ylim(1, 14)
ax.set_xlim(np.min(x), np.max(x))

ax.plot(x, y, "r-") # Plot of the data use lines

#add title, axis titles, and legend
ax.set_title("Acid titration")
ax.set_xlabel("Moles NaOH")
ax.set_ylabel("pH")
#ax.legend(("y data"), loc='upper left')

plt.show()


Based on the answer, here is what I came up with. Any other suggestions?

# my pH calculation class
# assume two distinct pKa's  solution is a quartic equation
class pH:
#things that don't change
ka1 = 2.479496e-6
ka2 = 1.87438e-9
kw = 1.019230e-14
kSum = ka1+ka2
kProd = ka1*ka2
e = -kw*kProd

#things that only depend on Facid1 and Facid2
def __init__(self, Facid1, Facid2):
self.c = -(self.ka1*Facid1+self.ka2*Facid2)+self.kProd-self.kw
self.d = self.kProd*(Facid1+Facid2)+self.kw*(self.kSum)

#only calculate things that depend on base
def pHCalc(self, base):
pMatrix = [[0, 0, 0, -self.e],  #construct the companion matrix
[1, 0, 0, self.d-base*self.kProd],
[0, 1, 0, -(self.c+self.kSum*base)],
[0, 0, 1, -(self.kSum+base)]]
myVals = la.eigvals(pMatrix)
return -np.log10(np.max(myVals)) #need the one positive root


NumPy computes the roots of a polynomial by first constructing the companion matrix in Python and then solving the eigenvalues with LAPACK. The companion matrix case looks like this using your variables (as a==1):
[0 0 0 -e

You should be able to save some time by updating a matrix like this directly on each iteration of base. Then use numpy.linalg.eigvals(m).max() to obtain the largest eigenvalue. See the sources.