Python Advanced: Polynomial Class

Polynomials



Introduction

Beauty in Math by Drew Shannon

If you have been to highschool, you will have encountered the terms polynomial and polynomial function. This chapter of our Python tutorial is completely on polynomials, i.e. we will define a class to define polynomials. The following is an example of a polynomial with the degree 4:

$$p(x) = x^4 - 4 \cdot x^2 + 3 \cdot x$$

You will find out that there are lots of similarities to integers. We will define various arithmetic operations for polynomials in our class, like addition, subtraction, multiplication and division. Our polynomial class will also provide means to calculate the derivation and the integral of polynomials. We will not miss out on plotting polynomials.

There is a lot of beaty in polynomials and above all in how they can be implemented as a Python class. We like to say thanks to Drew Shanon who granted us permission to use his great picture, treating math as art!

Short Mathematical Introduction

We will only deal with polynomial in a single indeterminate (also called variable) x. A general form of a polynomial in a single indeterminate looks like this:

$$a_n \cdot x^n + a_{n-1} \cdot x^{n-1} + \ldots + a_2 \cdot x^2 + a_1 \cdot x + a_0$$

where $a_0, a_1, ... a_n$ are the constants - non-negative integers - and $x$ is the indeterminate or variable. The term "indeterminate" means that $x$ represents no particular value, but any value may be substituted for it.

This expression is usually written with the summation operator:

$$\sum_{k=0}^{n} a_k \cdot x^k = a_n \cdot x^n + a_{n-1} \cdot x^{n-1} + \ldots + a_2 \cdot x^2 + a_1 \cdot x + a_0$$

A polynomial function is a function that can be defined by evaluating a polynomial. A function f of one argument can be defined as:

$$f(x) = \sum_{k=0}^{n} a_k \cdot x^k$$

Polynomial Functions with Python

It's easy to implement polynomial functions in Python. As an example we define the polynomial function given in the introduction of this chapter, i.e. $p(x) = x^4 - 4 \cdot x^2 + 3 \cdot x$

The Python code for this polynomial function looks like this:

In [6]:
def p(x):
    return x**4 - 4*x**2 + 3*x

We can call this function like any other function:

In [7]:
for x in [-1, 0, 2, 3.4]:
    print(x, p(x))
-1 -6
0 0
2 6
3.4 97.59359999999998

In [8]:
import numpy as np
import matplotlib.pyplot as plt

X = np.linspace(-3, 3, 50, endpoint=True)
F = p(X)
plt.plot(X,F)

plt.show()

Polynomial Class

We will define now a class for polynomial functions. We will build on an idea which we have developed in the chapter on decorators of our Python tutorial. We introduced polynomial factories.

A polynomial is uniquely determined by its coefficients. This means, an instance of our polynomial class needs a list or tuple to define the coefficients.

In [9]:
class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        # for reasons of efficiency we save the coefficients in reverse order,
        # i.e. a_0, a_1, ... a_n
        self.coefficients = coefficients[::-1] # tuple is also turned into list
     
    # method to return the canonical string representation of a polynomial
    def __repr__(self):
        return "Polynomial" + str(self.coefficients)
        

We can the instantiate the polynomial of our previous example polynomial function like this:

In [10]:
p = Polynomial(4, 0, -4, 3, 0)
print(p)
Polynomial(0, 3, -4, 0, 4)

So far, we have defined polynomials, but what we actually need are polynomial functions. For this purpose, we turn instances of the Polynomial class into callables by defining the call method:

In [11]:
class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        # for reasons of efficiency we save the coefficients in reverse order,
        # i.e. a_0, a_1, ... a_n
        self.coefficients = coefficients[::-1] # tuple is also turned into list
     
    # method to return the canonical string representation of a polynomial
    def __repr__(self):
        return "Polynomial" + str(self.coefficients)
    
    def __call__(self, x):    
        res = 0
        for index, coeff in enumerate(self.coefficients):
            res += coeff * x** index
        return res 

It is possible now to call an instance with an argument. This means that an instance behaves like a polynomial function:

In [12]:
p = Polynomial(4, 0, -4, 3, 0)

for x in range(-3, 3):
    print(x, p(x))
-3 279
-2 42
-1 -3
0 0
1 3
2 54

In [13]:
import matplotlib.pyplot as plt

X = np.linspace(-3, 3, 50, endpoint=True)
F = p(X)
plt.plot(X,F)

plt.show()

It is possible to define addition and subtractions for polynomials. All we have to do is to add or subtract the coefficients with the same exponents from both polynomials.

If we have polynomial functions $f(x) = \sum_{k=0}^{n} a_k \cdot x^k$ and $g(x) = \sum_{k=0}^{n} b_k \cdot x^k$, the addition is defined as

$$(f+g)(x) = \sum_{k=0}^{n} (a_k + b_k) \cdot x^k$$

and correspondingly subtraction is defined as

$$(f-g)(x) = \sum_{k=0}^{n} (a_k - b_k) \cdot x^k$$

Before we can add the methods __add__ and __sub__, which are necessary for addition and subtraction, we add a generatur zip_longest(). It works similar to zip for two parameters, but it does not stop if one of the iterators is exhausted but uses "fillchar" instead.

In [14]:
def zip_longest(iter1, iter2, fillchar=None):
    
    for i in range(max(len(iter1), len(iter2))):
        if i >= len(iter1):
            yield (fillchar, iter2[i])
        elif i >= len(iter2):
            yield (iter1[i], fillchar)
        else:
            yield (iter1[i], iter2[i])
        i += 1

p1 = (2,)
p2 = (-1, 4, 5)
for x in zip_longest(p1, p2, fillchar=0):
    print(x)
(2, -1)
(0, 4)
(0, 5)

We will add this generator to our class as a static method. We are capable of adding the __add__ and __sub__ methods now as well:

In [15]:
import numpy as np
import matplotlib.pyplot as plt

class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        # for reasons of efficiency we save the coefficients in reverse order,
        # i.e. a_0, a_1, ... a_n
        self.coefficients = coefficients[::-1] # tuple is also turned into list
     
    # method to return the canonical string representation of a polynomial
    def __repr__(self):
        return "Polynomial" + str(self.coefficients)
    
    def __call__(self, x):    
        res = 0
        for index, coeff in enumerate(self.coefficients):
            res += coeff * x** index
        return res 
    
    def degree(self):
        return len(self.coefficients)
    
    @staticmethod
    def zip_longest(iter1, iter2, fillchar=None):    
        for i in range(max(len(iter1), len(iter2))):
            if i >= len(iter1):
                yield (fillchar, iter2[i])
            elif i >= len(iter2):
                yield (iter1[i], fillchar)
            else:
                yield (iter1[i], iter2[i])
            i += 1   
            
    def __add__(self, other):
        c1 = self.coefficients
        c2 = other.coefficients
        res = [sum(t) for t in Polynomial.zip_longest(c1, c2)]
        return Polynomial(*res)
    
    def __sub__(self, other):
        c1 = self.coefficients
        c2 = other.coefficients
        
        res = [t1-t2 for t1, t2 in Polynomial.zip_longest(c1, c2)]
        return Polynomial(*res)
    

            
p1 = Polynomial(4, 0, -4, 3, 0)
p2 = Polynomial(-0.8, 2.3, 0.5, 1, 0.2)

p_sum = p1 + p2
p_diff = p1 - p2

X = np.linspace(-3, 3, 50, endpoint=True)
F1 = p1(X)
F2 = p2(X)
F_sum = p_sum(X)
F_diff = p_diff(X)
plt.plot(X, F1, label="F1")
plt.plot(X, F2, label="F2")
plt.plot(X, F_sum, label="F_sum")
plt.plot(X, F_diff, label="F_diff")

plt.legend()
plt.show()

It is incredibly easy to add differentiation to our class. Mathematically, it is defined as

$$f'(x) = \sum_{k=0}^{n} k \cdot a_k \cdot x^{k-1}$$

if
$$f(x) = \sum_{k=0}^{n} a_k \cdot x^k$$

This can be easily implemented in our method 'derivative':

In [78]:
import numpy as np
import matplotlib.pyplot as plt

class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        # for reasons of efficiency we save the coefficients in reverse order,
        # i.e. a_0, a_1, ... a_n
        self.coefficients = coefficients[::-1] # tuple is also turned into list
     
    # method to return the canonical string representation of a polynomial
    def __repr__(self):
        return "Polynomial" + str(self.coefficients)
    
    def __call__(self, x):    
        res = 0
        for index, coeff in enumerate(self.coefficients):
            res += coeff * x** index
        return res 
    
    def degree(self):
        return len(self.coefficients)
    
    @staticmethod
    def zip_longest(iter1, iter2, fillchar=None):    
        for i in range(max(len(iter1), len(iter2))):
            if i >= len(iter1):
                yield (fillchar, iter2[i])
            elif i >= len(iter2):
                yield (iter1[i], fillchar)
            else:
                yield (iter1[i], iter2[i])
            i += 1   
            
    def __add__(self, other):
        c1 = self.coefficients
        c2 = other.coefficients
        res = [sum(t) for t in Polynomial.zip_longest(c1, c2)]
        return Polynomial(*res)
    
    def __sub__(self, other):
        c1 = self.coefficients
        c2 = other.coefficients
        
        res = [t1-t2 for t1, t2 in Polynomial.zip_longest(c1, c2)]
        return Polynomial(*res)
    
    def derivative(self):
        derived_coeffs = []
        exponent = 1
        for i in range(1, len(self.coefficients)):
            derived_coeffs.append(self.coefficients[i] * exponent)
            exponent += 1
        return Polynomial(*derived_coeffs)
    
    def __str__(self):
        res = ""
        for i in range(len(self.coefficients)-1, -1, -1):
            res +=  str(self.coefficients[i]) + "x^" + str(i) + " + "
        if res.endswith(" + "):
            res = res[:-3]
        return res

            
p = Polynomial(-0.8, 2.3, 0.5, 1, 0.2)

p_der = p.derivative()

X = np.linspace(-2, 3, 50, endpoint=True)

F = p(X)
F_derivative = p_der(X)
plt.plot(X, F, label="F")
plt.plot(X, F_derivative, label="F_der")

plt.legend()
plt.show()
In [ ]: