Some Links

Tutorials

License

Computations with Polynomials

Numpy is useful anywhere complex mathematical computations are
involved. This includes engineering problems for which Matlab, a
very expensive product, is often used. It is easy to accomplish all
that work using Numpy as well. See
http://www.scipy.org/NumPy_for_Matlab_Users
for more details. In this article, we will explore some of the common
computational needs and how they can be met using Numpy.

Regression Analysis/Least Square Fit

A very common requirement you may have is that you have collected
a lot of data and have a model. You need to test the whether the data
fits the model. In this article, you can explore the Galilean model
because as mentioned in
http://www.vpri.org/pdf/rn2005001_learning.pdf
:“This subject has been extensively studied with college students
in US: 70% (including science majors) fail to understand this
Galilean model of gravity near the surface of the Earth.” It
roughly means that a large proportion of the population believes that
heavier objects fall faster. With the technology available these
days, experimental testing of these beliefs is not difficult.

You can take a video of a few falling objects and convert them
into frames. So, you may get about 30 frames per second. Identifying
an object in a video frame will not be very precise but the errors in
measurement are an inherent part of research. The model you wish to
verify is

y
= ½ g t2 + v0 t + y0

You have a list of positions of the object every 1/30th
of a second. You can fit the values to a 2nd degree
polynomial least square fit algorithm. Here is the code you would
use:

import numpy as np

values =
[195.7,193.51,190.59,185.84,180.73,173.43,

165.76,157.36,147.87,136.92,124.87]

y = np.array(values)

t = np.linspace(0, 10.0/30, 11)

poly_coeff = np.polyfit(t, y, 2)

print poly_coeff

print 'g = ', 2*poly_coeff[0]

You convert a Python list of 11 values into a numpy array. The
linspace method creates an equally spaced array with first value
being 0 and the last (11th) value being 1/3 sec. The
polyfit method will fit the data to a polynomial with the degree
being 2 and return an array of 3 numbers. These three numbers are the
coefficients of the polynomial starting with the highest power first.
So, the first value is half of g. The above data gives g as -963
cm/sec2 . The negative sign indicating downward direction.

How good was the fit?

Before discussing that, let us look at the relationship between an
array and a polynomial.

Polynomials and Arrays

Any polynomial can be represented by a vector by just keeping
track of the coefficients. E.g.

[1,0,-4] would represent (x2 – 4). Try the following
examples:

>>> import numpy as np

>>> p2 =
np.poly1d([1,0,-4])

>>> p2(0)

>>> p2([0,1,2])

You have created a object of the type polynomial (or simply a
polynomial function) of degree 2. Now, you can use it like you would
a function. The result p2(0) would be -4 and for p2([0,1,2]) would be
the array, [-4, -3, 0].

The best way to see how well the polynomial you obtained fits the
data is to visually see it through a plot.

Matplotlib

The python-matplotlib package must be installed.

You will want to plot the original data and the values computed
from the polynomial. Hence, add the following code after obtaining
the least square fit to your data:

import matplotlib.pyplot as plt

plt.plot(t, y, 'o')

plt.plot(t,
np.poly1d(poly_coeff)(t), '-')

plt.show()

The first call to the plot function plots small circles at the
actual data points pairs (ti,yi).

In the second call to the plot function, you are creating a
polynomial function from the array poly_coeff and, then computing the
values of this function for each value of t. Finally, you are
plotting these values against t as a line. The resulting fit for the
above example is pretty good as can be seen in the following figure:

The Roots

Solving a polynomial is another common requirement. The roots of
the polynomial fit to the data would indicate when the object would
be on the ground. Add the following lines to your code:

roots = np.roots(poly_coeff)

print 'Roots', roots

The output of the program, including the earlier code, would be:

[-481.72027972 -52.31748252
195.86951049]

g = -963.440559441

Roots [-0.69426607 0.58566055]

The meaningful root in this case is .585 sec. It may not be very
useful in this context but it does become valuable when dropping food
packets or bombs from an aeroplane.

A polynomial can be represented by the coefficients as we have
discussed. It can equally well be represented by its roots. So, if
you create an object of type polynomial, both representations are
easily available to you, as follows

>>>
pf=np.poly1d(poly_coeff)

>>> print 'The roots',
pf.r

The roots [-0.69426607
0.58566055]

>>> print 'The
coefficients', pf.c

The coefficients [-481.72027972
-52.31748252 195.86951049]

Multiple Data Sets

Typically, there will be multiple runs and you will need to
analyse each set of data. In the following example, you have two runs
of data – one for a stone and second for a small plastic cap. Try
the following code:

import numpy as np

set1 =
[195.7,193.51,190.59,185.84,180.73,173.43,

165.76,157.36,147.87,136.92,124.87]

set2 =
[194.04,192.6,189.36,185.04,180.36,173.52,

167.76,158.76,148.68,138.6,127.8]

y = np.array([set1,set2])

print 'Shape of the array',
y.shape

t = np.linspace(0, 10.0/30, 11)

res = np.polyfit(t,
y.transpose(), 2)

print 'The shape of the result',
res.shape

print 'g = ', 2*res[0]

The array y is list of data values for each set. However, polyfit
requires list to grouped by corresponding values for all trial data
sets. So, you need to take the transpose of the array y to fit the
needs of the polyfit funtion. The result will be a set of 3
coefficients, one for each set of data. You get a pair of values for
g. The output will now be:

Shape of the array (2, 11)

The shape of the result (3, 2)

g = [-963.44055944
-953.11888112]

You may need to do more experiments with plastic caps! But
meanwhile, you can plot the two sets of data and the fitted lines as
below:

import matplotlib.pyplot as plt

plt.plot(t,set1,'o',
label='Dataset 1')

plt.plot(t,set2,'d',
label='Dataset 2')

coeff = res.transpose()

plt.plot(t,np.poly1d(coeff[0])(t),'-',
label='Fit 1')

plt.plot(t,np.poly1d(coeff[1])(t),'-',
label='Fit 2')

plt.legend()

plt.show()

Pyplot will choose a different color for each plot. The first set
of points are drawn as small circles. The second set of points are
drawn as little diamonds. The last two are continuous lines. Note
that you needed to take the transpose of the result matrix for easy
access to the coefficients of each polynomial. A legend will be
useful in this case. See the figure below.

You can find the details of all the
options for plotting by using inbuilt help of python:

>>> help(plt.plot)

What Else

Polynomials are very useful. The integral of a polynomial is
another polynomial. The derivative is another polynomial. Numpy
includes functions, polyint and polyder for this purpose. The
function polymul will return another polynomial which is the product
of two polynomials. There are a wealth of functions available in
numpy which makes working with polynomials in education and research
a much easier task.

Scipy module of Python includes everything of numpy and some more.
So, you can replace the import statement by 'import scipy as np' and
your code will still run fine. In case you are still hesitating to
use Python for Science, the article “Python in Science: How long
until a Nobel Prize? | And Now For Something Completely
Different” by Doug Hellmann might interest you -
http://www.doughellmann.com/articles/CompletelyDifferent-2007-11-science/