### 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.

(A fun aside: Hammer and Feather Drop - Apollo 15 astronaut David Scott on the Moon demonstrates that Galileo was right. http://en.wikipedia.org/wiki/Galileo%27s_Leaning_Tower_of_Pisa_experiment )

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 t````2 ``+ v``0`` t + y``0 `

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`

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`

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)(t),'-', label='Fit 1')```

```plt.plot(t,np.poly1d(coeff)(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/