Polynome / Polynomklasse



Einführung

Beauty in Math by Drew Shannon

Viele haben bereits in der Schule die Begriffe Polynom und Polynomfunktion kennengelernt. In diesem Kapitel unseres Python-Kurses geht es um Polynome, d.h. wir werden eine Klasse definieren, um Polynome zu definieren. Nachfolgend ein Beispiel für ein Polynom mit dem Grad 4:

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

Wir werden sehen, dass man mit Polynomen so rechnen kann wie mit Integers und Floats auch. Operationen, die auch in Python implementiert werden können. In diesem Kurs werden wir verschiedene arithmetische Operationen für Polynome definieren, wie Addition, Subtraktion, Multiplikation und Division. Unsere Polynomklasse wird auch Mittel zur Berechnung der Ableitung und des Integrals von Polynomen bereitstellen. Auch das Plotten von Polynomen wird nicht zu kurz kommen.

Es gibt eine Menge Anmut oder Schönheit in Polynomen und vor allem darin, wie sie als Python-Klasse implementiert werden können. Wir bedanken uns bei Drew Shanon, der uns die Erlaubnis erteilt hat, sein großartiges Bild zu verwenden, das Mathe als Kunst behandelt!

Kurze mathematische Einführung

Wir werden uns nur mit Polynomen in einer einzigen Unbestimmten (Variablen) x beschäftigen. Eine allgemeine Form eines Polynoms in einer einzigen unbestimmten Zahl sieht so aus:

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

wobei $a_0, a_1, ... a_n$ die Konstanten des Polynoms sind und $x$ die Unbestimmte oder Variable ist. Der Begriff "unbestimmt" bedeutet, dass $x$ keinen bestimmten Wert darstellt, sondern durch einen beliebigen Wert ersetzt werden kann.

Dieser Ausdruck wird gewöhnlich mit dem Summationsoperator geschrieben:

$$\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$$

Eine Polynomfunktion ist eine Funktion, die durch Auswertung eines Polynoms definiert werden kann. Eine Funktion f mit einem Argument kann definiert werden als:

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

Polynom-Funktionen mit Python

Es ist einfach, Polynomfunktionen in Python zu implementieren. Als Beispiel definieren wir die Polynomfunktion aus der Einleitung dieses Kapitels, d.h. $p(x) = x^4 - 4 \cdot x^2 + 3 \cdot x$

Der Python-Code für diese Polynomfunktion sieht folgendermaßen aus:

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

Wir können diese Funktion wie jede andere Funktion aufrufen:

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

Plotten können wir die Funktion auf einfache Art und Weise mit Matplotlib:

In [3]:
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()
No description has been provided for this image

Polynom-Klasse in Python

Wir werden nun eine Klasse für Polynom-Funktionen in Python definieren. Dabei bauen wir auf einer Idee auf, die wir im Kapitel über Dekorateure in unserem Python-Tutorial entwickelt haben.

Ein Polynom ist eindeutig durch seine Koeffizienten bestimmt. Das heißt, eine Instanz unserer Polynomklasse benötigt eine Liste oder ein Tupel, um die Koeffizienten zu definieren.

In [4]:
class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list
     
    def __repr__(self):
        """
        method to return the canonical string representation 
        of a polynomial.
        """
        return "Polynomial" + str(tuple(self.coefficients))
        

Wir können das Polynom unserer vorherigen Beispiel-Polynomfunktion wie folgt instanziieren:

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

p2 = eval(repr(p))
print(p2)
Polynomial(1, 0, -4, 3, 0)
Polynomial(1, 0, -4, 3, 0)

Nachdem wir die kanonische String-Darstellung eines Polynoms mit der repr-Funktion definiert haben, wollen wir auch eine Ausgabe-Methode definieren, die verwendet werden kann, um eine für Menschen freundlichere Version zu erstellen. Wir schreiben eine LaTex-Darstellung, die verwendet werden kann, um die Funktion auf eine schöne Weise auszugeben:

In [6]:
class Polynomial:
    
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list
     
    
    def __repr__(self):
        """
        method to return the canonical string representation 
        of a polynomial.
        """
        return "Polynomial" + str(tuple(self.coefficients))

    
    def __str__(self):
        
        def x_expr(degree):
            if degree == 0:
                res = ""
            elif degree == 1:
                res = "x"
            else:
                res = "x^"+str(degree)
            return res

        degree = len(self.coefficients) - 1
        res = ""

        for i in range(0, degree+1):
            coeff = self.coefficients[i]
            # nothing has to be done if coeff is 0:
            if abs(coeff) == 1 and i < degree:
                # 1 in front of x shouldn't occur, e.g. x instead of 1x
                # but we need the plus or minus sign:
                res += f"{'+' if coeff>0 else '-'}{x_expr(degree-i)}"  
            elif coeff != 0:
                res += f"{coeff:+g}{x_expr(degree-i)}" 

        return res.lstrip('+')    # removing leading '+'
In [7]:
polys = [Polynomial(1, 0, -4, 3, 0),
         Polynomial(2, 0),
         Polynomial(4, 1, -1),
         Polynomial(3, 0, -5, 2, 7),
         Polynomial(-42)]

# Ausgabe geeignet zur Verwendung in LaTeX:
for count, poly in enumerate(polys):
    print(f"$p_{count} = {str(poly)}$")
$p_0 = x^4-4x^2+3x$
$p_1 = 2x$
$p_2 = 4x^2+x-1$
$p_3 = 3x^4-5x^2+2x+7$
$p_4 = -42$

Benutzen wir dies in LaTeX, sieht die Ausgabe wie folgt aus:

$p_0 = x^4-4x^2+3x$

$p_1 = 2x$

$p_2 = 4x^2+x-1$

$p_3 = 3x^4-5x^2+2x+7$

$p_4 = -42$

Bislang haben wir Polynome definiert, aber was wir eigentlich brauchen, sind Polynomfunktionen. Zu diesem Zweck verwandeln wir Instanzen der Klasse Polynomial in Callabes, also aufrufbare Klasseninstanzen. Dazu definieren wir Methode call:

In [8]:
class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list
    
    # The __repr__ and __str__ method can be included here,
    # but is not necessary for the immediately following code
            
    def __call__(self, x):    
        res = 0
        for index, coeff in enumerate(self.coefficients[::-1]):
            res += coeff * x** index
        return res 

Es ist nun möglich, eine Instanz unserer Klasse wie eine Funktion aufzurufen. Wir rufen sie mit einem Argument auf und die Instanz, die ein Callable ist, verhält sich wie eine Polynomfunktion:

In [9]:
p = Polynomial(3, 0, -5, 2, 1)
print(p)

for x in range(-3, 3):
    print(x, p(x))
<__main__.Polynomial object at 0x7fc61d9d4ee0>
-3 193
-2 25
-1 -3
0 1
1 1
2 33

Nun wollen wir die vorher definierte Funkton ausgeben:

In [10]:
import matplotlib.pyplot as plt

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

plt.show()
No description has been provided for this image

Bevor wir unsere Klasse weiter verfeinern, wollen wir eine numerisch effizientere Variante für die Berechnung von Polynomen verwenden. Wir werden diesen Algorithmus in unserer Methode __call__ verwenden. Wir haben diese Variante in unserem Kapitel über Dekorateure vorgestellt.

Jedes Polynom

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

kann auch in folgender Form geschrieben werden

$f(x) = (a_n x + x_{n-1})x + \cdot + a_1)x + a_0$

Um empirisch zu sehen, dass sie äquivalent sind, schreiben wir unsere Klassendefinition um, nennen sie aber Polynom2, so dass wir beide Versionen verwenden können:

In [11]:
class Polynomial2:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list
         
    # The __repr__ and __str__ method can be included here,
    # but is not necessary for the immediately following code
    
    def __call__(self, x):    
        res = 0
        for coeff in self.coefficients:
            res = res * x + coeff
        return res 

Wir testen beide Klassen, indem wir die Ergebnisse vergleichen:

In [12]:
p1 = Polynomial(-4, 3, 0)
p2 = Polynomial2(-4, 3, 0)
res = all((p1(x)==p2(x) for x in range(-10, 10)))
print(res)
True

Es ist möglich, Additionen und Subtraktionen für Polynome zu definieren. Alles, was wir tun müssen, ist, die Koeffizienten mit den gleichen Exponenten von beiden Polynomen zu addieren oder zu subtrahieren.

Wenn wir Polynomfunktionen $f(x) = \sum_{k=0}^{n} a_k \cdot x^k$ und $g(x) = \sum_{k=0}^{m} b_k \cdot x^k$ haben, ist die Addition definiert als

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

und dementsprechend ist die Subtraktion definiert als

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

Bevor wir die Methoden __add__ und __sub__, die für Addition und Subtraktion notwendig sind, hinzufügen können, fügen wir einen Generator zip_longest() hinzu. Sie funktioniert ähnlich wie zip für zwei Parameter, hält aber nicht an, wenn einer der Iteratoren erschöpft ist, sondern verwendet stattdessen "fillvalue".

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

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

Den eben definierten Generator zip_longest benötigen wir jedoch nicht, da das Modul itertools einen Generator zip_longest enthält.

Jetzt sind wir in der Lage, auch die Methoden __add__ und __sub__ unserer Klasse hinzuzufügen.

zip_longest nimmt eine beliebige Anzahl von iterierbaren Objekten und einen Schlüsselwort-Parameter fillvalue. Es erzeugt einen Iterator, der die Elemente der übergebenen iterierbaren Objekte zusammenfasst. Wenn die iterierbaren Objekte nicht die gleiche Länge haben, werden fehlende Werte mit dem Wert von fillvalue aufgefüllt. Die Iteration wird fortgesetzt, bis die längste Iterabel erschöpft ist.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from itertools import zip_longest

class Polynomial:
    
    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list
     
    def __repr__(self):
        """
        method to return the canonical string representation 
        of a polynomial.
   
        """
        return "Polynomial" + eval(str_s)
            
    def __call__(self, x):    
        res = 0
        for coeff in self.coefficients:
            res = res * x + coeff
        return res 
    
    def degree(self):
        return len(self.coefficients) - 1
            
    def __add__(self, other):
        c1 = self.coefficients[::-1]
        c2 = other.coefficients[::-1]
        res = [sum(t) for t in zip_longest(c1, c2, fillvalue=0)]
        return Polynomial(*res[::-1])
    
    def __sub__(self, other):
        c1 = self.coefficients[::-1]
        c2 = other.coefficients[::-1]
        
        res = [t1-t2 for t1, t2 in zip_longest(c1, c2, fillvalue=0)]
        return Polynomial(*res[::-1])
    

            
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()
No description has been provided for this image

Es ist unglaublich einfach, auch die Berechnung der Ableitung in unserer Klassendefinition mit aufzunehmen. Mathematisch ist sie definiert als

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

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

Dies kann leicht in unserer Methode 'derivative' implementiert werden:

In [6]:
import numpy as np
import matplotlib.pyplot as plt
from itertools import zip_longest

class Polynomial:
 

    def __init__(self, *coefficients):
        """ input: coefficients are in the form a_n, ...a_1, a_0 
        """
        self.coefficients = list(coefficients) # tuple is turned into a list

        
    def __repr__(self):
        """
        method to return the canonical string representation 
        of a polynomial.
   
        """
        return "Polynomial" + str(tuple(self.coefficients))

    
    def __call__(self, x):    
        res = 0
        for coeff in self.coefficients:
            res = res * x + coeff
        return res 

    
    def degree(self):
        return len(self.coefficients)  - 1  

    
    def __add__(self, other):
        c1 = self.coefficients[::-1]
        c2 = other.coefficients[::-1]
        res = [sum(t) for t in zip_longest(c1, c2, fillvalue=0)]
        return Polynomial(*res[::-1])

    
    def __sub__(self, other):
        c1 = self.coefficients[::-1]
        c2 = other.coefficients[::-1]
        
        res = [t1-t2 for t1, t2 in zip_longest(c1, c2, fillvalue=0)]
        return Polynomial(*res)
 

    def derivative(self):
        derived_coeffs = []
        exponent = len(self.coefficients) - 1
        for i in range(len(self.coefficients)-1):
            derived_coeffs.append(self.coefficients[i] * exponent)
            exponent -= 1
        return Polynomial(*derived_coeffs)

    
    def __str__(self):
        
        def x_expr(degree):
            if degree == 0:
                res = ""
            elif degree == 1:
                res = "x"
            else:
                res = "x^"+str(degree)
            return res

        degree = len(self.coefficients) - 1
        res = ""

        for i in range(0, degree+1):
            coeff = self.coefficients[i]
            # nothing has to be done if coeff is 0:
            if abs(coeff) == 1 and i < degree:
                # 1 in front of x shouldn't occur, e.g. x instead of 1x
                # but we need the plus or minus sign:
                res += f"{'+' if coeff>0 else '-'}{x_expr(degree-i)}"  
            elif coeff != 0:
                res += f"{coeff:+g}{x_expr(degree-i)}" 

        return res.lstrip('+')    # removing leading '+'

            
In [7]:
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()
No description has been provided for this image

Nun benutzen wir die Addition und Subtraktion von Polynomen:

In [8]:
p = Polynomial(1, 2, -3, 4, -5)

p2 = Polynomial(3, 2, -3, 4, 3)

(p + p2)
Out[8]:
Polynomial(4, 4, -6, 8, -2)
In [9]:
import matplotlib.pyplot as plt

plt.plot(X, p(X), label="p")
plt.plot(X, p2(X), label="p2")
plt.plot(X, (p + p2)(X), label="p + p2")
plt.plot(X, (p - p2)(X), label="p - p2")

plt.legend()
plt.show()
No description has been provided for this image