p(x) = 1 + 2x + 3x^2

p(x) = 1 T_0(x) + 2 T_1(x) + 3 T_2(x)

p(x) = \sum_{i=0}^n c_i T_i(x)

Basics¶

First we need a polynomial class and a polynomial instance to play with.The classes can be imported directly from the polynomial package or fromthe module of the relevant type. Here we import from the package and usethe conventional Polynomial class because of its familiarity:

>>>fromnumpy.polynomialimportPolynomialasP>>>p=P([1,2,3])>>>pPolynomial([ 1.,  2.,  3.], domain=[-1,  1], window=[-1,  1])

Note that there are three parts to the long version of the printout. Thefirst is the coefficients, the second is the domain, and the third is thewindow:

>>>p.coefarray([ 1.,  2.,  3.])>>>p.domainarray([-1.,  1.])>>>p.windowarray([-1.,  1.])

Printing a polynomial yields a shorter form without the domainand window:

>>>printppoly([ 1.  2.  3.])

We will deal with the domain and window when we get to fitting, for the momentwe ignore them and run through the basic algebraic and arithmetic operations.

Addition and Subtraction:

>>>p+pPolynomial([ 2.,  4.,  6.], domain=[-1,  1], window=[-1,  1])>>>p-pPolynomial([ 0.], domain=[-1,  1], window=[-1,  1])

Multiplication:

>>>p*pPolynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])

Powers:

>>>p**2Polynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])

Division:

Floor division, ‘//’, is the division operator for the polynomial classes,polynomials are treated like integers in this regard. For Python versions <3.x the ‘/’ operator maps to ‘//’, as it does for Python, for laterversions the ‘/’ will only work for division by scalars. At some point itwill be deprecated:

>>>p//P([-1,1])Polynomial([ 5.,  3.], domain=[-1,  1], window=[-1,  1])

Remainder:

>>>p%P([-1,1])Polynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Divmod:

>>>quo,rem=divmod(p,P([-1,1]))>>>quoPolynomial([ 5.,  3.], domain=[-1,  1], window=[-1,  1])>>>remPolynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Evaluation:

>>>x=np.arange(5)>>>p(x)array([  1.,   6.,  17.,  34.,  57.])>>>x=np.arange(6).reshape(3,2)>>>p(x)array([[  1.,   6.],       [ 17.,  34.],       [ 57.,  86.]])

Substitution:

Substitute a polynomial for x and expand the result. Here we substitutep in itself leading to a new polynomial of degree 4 after expansion. Ifthe polynomials are regarded as functions this is composition offunctions:

>>>p(p)Polynomial([  6.,  16.,  36.,  36.,  27.], domain=[-1,  1], window=[-1,  1])

Roots:

>>>p.roots()array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])

It isn’t always convenient to explicitly use Polynomial instances, sotuples, lists, arrays, and scalars are automatically cast in the arithmeticoperations:

>>>p+[1,2,3]Polynomial([ 2.,  4.,  6.], domain=[-1,  1], window=[-1,  1])>>>[1,2,3]*pPolynomial([  1.,   4.,  10.,  12.,   9.], domain=[-1,  1], window=[-1,  1])>>>p/2Polynomial([ 0.5,  1. ,  1.5], domain=[-1,  1], window=[-1,  1])

Polynomials that differ in domain, window, or class can’t be mixed inarithmetic:

>>>fromnumpy.polynomialimportChebyshevasT>>>p+P([1],domain=[0,1])Traceback (most recent call last):  File"<stdin>", line1, in<module>  File"<string>", line213, in__add__TypeError:Domains differ>>>p+P([1],window=[0,1])Traceback (most recent call last):  File"<stdin>", line1, in<module>  File"<string>", line215, in__add__TypeError:Windows differ>>>p+T([1])Traceback (most recent call last):  File"<stdin>", line1, in<module>  File"<string>", line211, in__add__TypeError:Polynomial types differ

But different types can be used for substitution. In fact, this is howconversion of Polynomial classes among themselves is done for type, domain,and window casting:

>>>p(T([0,1]))Chebyshev([ 2.5,  2. ,  1.5], domain=[-1,  1], window=[-1,  1])

Which gives the polynomialp in Chebyshev form. This works becauseT_1(x) = x and substitutingx forx doesn’t changethe original polynomial. However, all the multiplications and divisionswill be done using Chebyshev series, hence the type of the result.

Calculus¶

Polynomial instances can be integrated and differentiated.:

>>>fromnumpy.polynomialimportPolynomialasP>>>p=P([2,6])>>>p.integ()Polynomial([ 0.,  2.,  3.], domain=[-1,  1], window=[-1,  1])>>>p.integ(2)Polynomial([ 0.,  0.,  1.,  1.], domain=[-1,  1], window=[-1,  1])

The first example integratesp once, the second example integrates ittwice. By default, the lower bound of the integration and the integrationconstant are 0, but both can be specified.:

>>>p.integ(lbnd=-1)Polynomial([-1.,  2.,  3.], domain=[-1,  1], window=[-1,  1])>>>p.integ(lbnd=-1,k=1)Polynomial([ 0.,  2.,  3.], domain=[-1,  1], window=[-1,  1])

In the first case the lower bound of the integration is set to -1 and theintegration constant is 0. In the second the constant of integration is setto 1 as well. Differentiation is simpler since the only option is thenumber of times the polynomial is differentiated:

>>>p=P([1,2,3])>>>p.deriv(1)Polynomial([ 2.,  6.], domain=[-1,  1], window=[-1,  1])>>>p.deriv(2)Polynomial([ 6.], domain=[-1,  1], window=[-1,  1])

Other Polynomial Constructors¶

Constructing polynomials by specifying coefficients is just one way ofobtaining a polynomial instance, they may also be created by specifyingtheir roots, by conversion from other polynomial types, and by leastsquares fits. Fitting is discussed in its own section, the other methodsare demonstrated below:

>>>fromnumpy.polynomialimportPolynomialasP>>>fromnumpy.polynomialimportChebyshevasT>>>p=P.fromroots([1,2,3])>>>pPolynomial([ -6.,  11.,  -6.,   1.], domain=[-1,  1], window=[-1,  1])>>>p.convert(kind=T)Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], domain=[-1,  1], window=[-1,  1])

The convert method can also convert domain and window:

>>>p.convert(kind=T,domain=[0,1])Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], [ 0.,  1.], [-1.,  1.])>>>p.convert(kind=P,domain=[0,1])Polynomial([-1.875,  2.875, -1.125,  0.125], [ 0.,  1.], [-1.,  1.])

In numpy versions >= 1.7.0 thebasis andcast class methods are alsoavailable. The cast method works like the convert method while the basismethod returns the basis polynomial of given degree:

>>>P.basis(3)Polynomial([ 0.,  0.,  0.,  1.], domain=[-1,  1], window=[-1,  1])>>>T.cast(p)Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], domain=[-1,  1], window=[-1,  1])

Conversions between types can be useful, but it isnot recommendedfor routine use. The loss of numerical precision in passing from aChebyshev series of degree 50 to a Polynomial series of the same degreecan make the results of numerical evaluation essentially random.

Fitting¶

Fitting is the reason that thedomain andwindow attributes are part ofthe convenience classes. To illustrate the problem, the values of the Chebyshevpolynomials up to degree 5 are plotted below.

>>>importmatplotlib.pyplotasplt>>>fromnumpy.polynomialimportChebyshevasT>>>x=np.linspace(-1,1,100)>>>foriinrange(6):ax=plt.plot(x,T.basis(i)(x),lw=2,label="$T_%d$"%i)...>>>plt.legend(loc="upper left")<matplotlib.legend.Legend object at 0x3b3ee10>>>>plt.show()

In the range -1 <=x <= 1 they are nice, equiripple functions lying between +/- 1.The same plots over the range -2 <=x <= 2 look very different:

>>>importmatplotlib.pyplotasplt>>>fromnumpy.polynomialimportChebyshevasT>>>x=np.linspace(-2,2,100)>>>foriinrange(6):ax=plt.plot(x,T.basis(i)(x),lw=2,label="$T_%d$"%i)...>>>plt.legend(loc="lower right")<matplotlib.legend.Legend object at 0x3b3ee10>>>>plt.show()

As can be seen, the “good” parts have shrunk to insignificance. In usingChebyshev polynomials for fitting we want to use the region wherex isbetween -1 and 1 and that is what thewindow specifies. However, it isunlikely that the data to be fit has all its data points in that interval,so we usedomain to specify the interval where the data points lie. Whenthe fit is done, the domain is first mapped to the window by a lineartransformation and the usual least squares fit is done using the mappeddata points. The window and domain of the fit are part of the returned seriesand are automatically used when computing values, derivatives, and such. Ifthey aren’t specified in the call the fitting routine will use the defaultwindow and the smallest domain that holds all the data points. This isillustrated below for a fit to a noisy sine curve.

>>>importnumpyasnp>>>importmatplotlib.pyplotasplt>>>fromnumpy.polynomialimportChebyshevasT>>>np.random.seed(11)>>>x=np.linspace(0,2*np.pi,20)>>>y=np.sin(x)+np.random.normal(scale=.1,size=x.shape)>>>p=T.fit(x,y,5)>>>plt.plot(x,y,'o')[<matplotlib.lines.Line2D object at 0x2136c10>]>>>xx,yy=p.linspace()>>>plt.plot(xx,yy,lw=2)[<matplotlib.lines.Line2D object at 0x1cf2890>]>>>p.domainarray([ 0.        ,  6.28318531])>>>p.windowarray([-1.,  1.])>>>plt.show()