问使用numpy多项式模块-有更好的方法吗？

EN

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

plt.show()

# 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