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()
ax = fig.add_subplot(111)

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

# Add my data
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.legend(("y data"), loc='upper left')


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

1 Answer 1


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
 1 0 0 -d
 0 1 0 -c
 0 0 1 -b]

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.

  • \$\begingroup\$ Yes, this works. I had already caught the need for a max(). I'm also thinking I could use a class construction to initialize some invariants, so I don't have to recalculate them each time. \$\endgroup\$ Commented Apr 3, 2013 at 20:22
  • 1
    \$\begingroup\$ Did some timing work with the full problem. Original code (using polynomial module) 5.48 sec, first step (using eigenvalue approach) 2.76 sec, final step (using eigenvalues and class as shown above) 2.68 sec. \$\endgroup\$ Commented Apr 3, 2013 at 21:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.