Performant Python

Polynomials

Overview:

  • Teaching: 5 min
  • Exercises: 15 min

Questions

  • What else can numpy do for me?

Objectives

  • Know that Numpy and SciPy provide a huge range of functionality
  • Know that Polynomial exist

Lectures on Numpy, Scipy

As well as a huge range of functionality there also exist lectures on use case for Numpy and the whole of Scipy. If you want to now want explore this material please do so. There is little value in trying to reproduce it here.

One example we will explore however is the polynomial as it gives an example of how numpy can help us produce very powerful efficient code efficiently. First what is numpy Polynomial module?

Polynomials

We can represent polynomials with the numpy class Polynomial from numpy.polynomial.Polynomial

  • Polynomial([a, b, c, d, e]) creates a polynomial $p(x) = a\,+\,b\,x \,+\,c\,x^2\,+\,d\,x^3\,+\,e\,x^4$.

For example:

  • Polynomial([1,2,3]) creates $p(x) = 1\,+\,2\,x \,+\,3\,x^2$

  • Polynomial([0,1,0,2,0,3]) creates $p(x) = x \,+\,2\,x^3\,+\,3\,x^5 $</p>

Let's import the class and have a quick play:

In [1]:
from numpy.polynomial import Polynomial as poly
In [2]:
my_poly = poly([1,2,3])
print(my_poly)
my_poly
poly([1. 2. 3.])
Out[2]:
$x \mapsto \text{1.0} + \text{2.0}\,x + \text{3.0}\,x^{2}$

Once created we can use help/autocompletion to explore/access the various functions available. We can access find the derivative of our polynomial:

In [3]:
my_poly.deriv()
Out[3]:
$x \mapsto \text{2.0} + \text{6.0}\,x$

If we assign this to a variable the the object is another polynomial:

In [4]:
my_deriv = my_poly.deriv()
print(my_deriv)
my_deriv
poly([2. 6.])
Out[4]:
$x \mapsto \text{2.0} + \text{6.0}\,x$
In [5]:
my_integ = my_poly.integ()
print(my_integ)
my_integ
poly([0. 1. 1. 1.])
Out[5]:
$x \mapsto \color{LightGray}{\text{0.0}} + \text{1.0}\,x + \text{1.0}\,x^{2} + \text{1.0}\,x^{3}$

Note that even numpy can't know the value of the constant when you integrate, though you can set it:

In [6]:
help(my_poly.integ)
Help on method integ in module numpy.polynomial._polybase:

integ(m=1, k=[], lbnd=None) method of numpy.polynomial.polynomial.Polynomial instance
    Integrate.
    
    Return a series instance that is the definite integral of the
    current series.
    
    Parameters
    ----------
    m : non-negative int
        The number of integrations to perform.
    k : array_like
        Integration constants. The first constant is applied to the
        first integration, the second to the second, and so on. The
        list of values must less than or equal to `m` in length and any
        missing values are set to zero.
    lbnd : Scalar
        The lower bound of the definite integral.
    
    Returns
    -------
    new_series : series
        A new series representing the integral. The domain is the same
        as the domain of the integrated series.

Calculting $\pi$ again

The Taylor series expansion for the trigonometric function $\arctan(y)$ is :

$$\arctan ( y) \, = \,y - \frac{y^3}{3} + \frac{y^5}{5} - \frac{y^7}{7} + \dots $$

Now, $\arctan(1) = \pi/4 $, so

$$ \pi = 4 \, \Big( 1- \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + ... \Big) $$

We can represent the series expansion using a numpy Polynomial, with coefficients:

$$0, +1, 0, -1/3, 0, +1/5, 0, -1/7,\ldots\,$$

and use it to approximate $\pi$.

Using this, write a code to estimate $\pi$. How does your code perform, how does it compare with the Monte Carlo approach?

Key Points:

  • Numpy and Scipy lectures are a valuable resource for learning to use their functionality
  • Making use of existing functionality will help you write correct, performant code more efficiently