18. Dez 2017
Lesedauer 10 Min.
Der Mathe-Turbo
Mathematik mit Python, Teil 2
Schneller rechnen mit numpy, scipy und sympy.
Im ersten Teil dieser zweiteiligen Serie zum Thema Mathematik mit Python [1] wurde die Array-Klasse im Modul numpy ausführlich vorgestellt. Diese Klasse beherbergt viele Array-Operationen, die das Programmieren erheblich vereinfachen. Dieser zweite Teil der Serie lenkt den Blick auf die mathematischen Lösungsverfahren im numpy-Modul: Statistik, Polynome, lineare Algebra und Matrizen-Rechnung. Außerdem
wird das Modul scipy vorgestellt, und auch etwas symbolische Mathematik mit sympy wird noch mit dabei sein.Bei den genannten mathematischen Verfahren ist eines immer ziemlich wichtig: die Performance. Wenn man ein Gleichungssystem mit einigen Hundert Unbekannten lösen muss, dann will man nicht stundenlang auf das Ergebnis warten. Da kann ich jedoch alle Leser beruhigen: Python ist sehr schnell, insbesondere beim Lösen von Rechenaufgaben! Das liegt daran, dass die Module numpy, scipy oder auch sympy in der Programmiersprache C implementiert sind. Dabei wurden sowohl die Belange von Rechnern mit begrenzten Systemressourcen (zum Beispiel Cache-Speicher) als auch erweiterte Hardware-Ressourcen (zum Beispiel Mehrkernprozessoren) berücksichtigt.
wird das Modul scipy vorgestellt, und auch etwas symbolische Mathematik mit sympy wird noch mit dabei sein.Bei den genannten mathematischen Verfahren ist eines immer ziemlich wichtig: die Performance. Wenn man ein Gleichungssystem mit einigen Hundert Unbekannten lösen muss, dann will man nicht stundenlang auf das Ergebnis warten. Da kann ich jedoch alle Leser beruhigen: Python ist sehr schnell, insbesondere beim Lösen von Rechenaufgaben! Das liegt daran, dass die Module numpy, scipy oder auch sympy in der Programmiersprache C implementiert sind. Dabei wurden sowohl die Belange von Rechnern mit begrenzten Systemressourcen (zum Beispiel Cache-Speicher) als auch erweiterte Hardware-Ressourcen (zum Beispiel Mehrkernprozessoren) berücksichtigt.
Statistik
Die statistischen Funktionen im Modul numpy sind ein gutes Einstiegsthema, da die Berechnungen überschaubar und gut nachvollziehbar sind. Mittelwert, Medianwert und gewichteter Durchschnitt sind leicht zu berechnen. Hier werden die FunkDabei wird ein zweidimensionales Array mithilfe des axis-Parameters in einer bestimmten Richtung verarbeitet und dann ein Ergebnisvektor zurückgegeben.Der Medianwert ist der nach dem Sortieren der Zahlen genau in der Mitte stehende Wert. Ist die Anzahl der Zahlen gerade, steht kein Wert genau in der Mitte, deshalb wird der Mittelwert zwischen den beiden mittleren Zahlen benutzt. Beim gewichteten Mittelwert muss der Funktion average ein Array mit Gewichtungen übergeben werden. Dieses Array hat die gleiche Größe wie das Daten-Array.Führen Sie das Programm aus Listing 1 aus, so erhalten Sie folgende Ergebnisse:Listing 1: Einfache Statistik
import numpy as <span class="hljs-built_in">np</span> <br/><br/>x = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">array</span>([<span class="hljs-number">1.0</span>, <span class="hljs-number">4.0</span>, <span class="hljs-number">9.0</span>, <span class="hljs-number">16.0</span>, <span class="hljs-number">25.0</span>]) <br/><br/># Normaler Mittelwert <br/>m1 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">mean</span>(x) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Mittelwert:"</span>, m1) <br/><br/># Medianwert <br/>m2 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">median</span>(x) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Medianwert:"</span>, m2) <br/><br/># Gewichteter Mittelwert <br/>m3 = <span class="hljs-built_in">np</span>.average(x, <br/> weights=[<span class="hljs-number">0.5</span>, <span class="hljs-number">1.0</span>, <span class="hljs-number">2.0</span>, <span class="hljs-number">1.0</span>, <span class="hljs-number">0.5</span>]) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Gewichtet: "</span>, m3)
<span class="hljs-selector-tag">Mittelwert</span>: 11<span class="hljs-selector-class">.0</span>
<span class="hljs-selector-tag">Medianwert</span>: 9<span class="hljs-selector-class">.0</span>
<span class="hljs-selector-tag">Gewichtet</span>: 10<span class="hljs-selector-class">.2</span>
Im nächsten Schritt sollen Standardabweichung und Varianz berechnet werden. Die Funktionen dafür heißen std und var. In Listing 2 wird ein Array mit Messwerten x definiert. Die Messwerte streuen um einen mittleren Wert. Die Standardabweichung gibt nun an, wie dicht die tatsächlichen Messwerte am mittleren Wert liegen. Die Varianz ist ebenfalls ein Maß für die Streuung, wird aber etwas anders berechnet.
Listing 2: Standardabweichung und Varianz
import numpy as np <br/><br/>x = np.array([<span class="hljs-number">5.2</span>, <span class="hljs-number">5.5</span>, <span class="hljs-number">4.8</span>, <span class="hljs-number">3.9</span>, <span class="hljs-number">4.5</span>, <span class="hljs-number">5.0</span>,<br/> <span class="hljs-number">5.2</span>, <span class="hljs-number">5.1</span>, <span class="hljs-number">4.8</span>, <span class="hljs-number">4.6</span>]) <br/><br/># <span class="hljs-number">1.</span> Mittelwert <br/>m = np.mean(x); <br/>print (<span class="hljs-string">"Mittelwert:"</span> , m) <br/><br/># <span class="hljs-number">2.</span> Standardabweichung <br/>s = np.std(x) <br/>print (<span class="hljs-string">"Standardabweichung:"</span>, s) <br/><br/># <span class="hljs-number">3.</span> Varianz <br/>v = np.var(x) <br/>print (<span class="hljs-string">"Varianz:"</span>, v) Polynome
Genug Statistik, jetzt sind die Polynome an der Reihe. Ein Polynom ist eine Potenzserie mit Koeffizienten ci, die das Polynom bestimmen:
P(<span class="hljs-keyword">x</span>)=<span class="hljs-keyword">c</span><sub><span class="hljs-number">0</span></sub> + <span class="hljs-keyword">c</span><sub><span class="hljs-number">1</span></sub><span class="hljs-keyword">x</span> + <span class="hljs-keyword">c</span><sub><span class="hljs-number">2</span></sub><span class="hljs-keyword">x</span><sup><span class="hljs-number">2</span></sup> + <span class="hljs-keyword">c</span><sub><span class="hljs-number">3</span></sub><span class="hljs-keyword">x</span><sup><span class="hljs-number">3</span></sup> +...+ c<sub>n</sub>x<sup>n</sup>
oder etwas einfacher hingeschrieben:
<span class="hljs-function"><span class="hljs-title">P</span><span class="hljs-params">(x)</span></span> = ∑ c<sub>i</sub>x<sup>i</sup>
Ein Beispiel für ein Polynom zweiten Grades wäre also:
P(<span class="hljs-keyword">x</span>) = <span class="hljs-number">5</span> – <span class="hljs-number">3</span><span class="hljs-keyword">x</span> + <span class="hljs-number">2</span><span class="hljs-keyword">x</span><sup><span class="hljs-number">2</span></sup>
Dieses Polynom kann man in Python ganz einfach definieren, indem man ein Polynomial-Objekt anlegt und die Koeffizienten als Array übergibt. Die Reihenfolge der Koeffizienten ist dabei zu beachten: [c0, c1, c2, c3, ...].
<span class="hljs-keyword">import</span> numpy <span class="hljs-keyword">as</span> np
from numpy.polynomial <span class="hljs-keyword">import</span> Polynomial
p = Polynomial([5.0, -3.0, 2.0])
print (p(4.5))
print (p(-1.0))
Nach der Definition des Polynomial-Objekts kann das Polynom P(x) mit beliebigen x-Werten berechnet werden. Das funktioniert natürlich auch direkt mit einem Array von x-Werten als Eingabe:
x = np.linspace(<span class="hljs-number">-5</span>, <span class="hljs-number">5</span>, <span class="hljs-number">11</span>)
y = p(x)
print (y)
Diese Zeilen erzeugen ein Array x mit elf Werten. Als Ergebnis erhalten Sie ein Array y mit folgenden elf Polynomwerten:
[ <span class="hljs-number">70.</span> <span class="hljs-number">49.</span> <span class="hljs-number">32.</span> <span class="hljs-number">19.</span> <span class="hljs-number">10.</span> <span class="hljs-number">5.</span> <span class="hljs-number">4.</span> <span class="hljs-number">7.</span> <span class="hljs-number">14.</span> <span class="hljs-number">25.</span> <span class="hljs-number">40.</span>]
Mit solchen Polynomial-Objekten können Sie nun auch algebraische Operationen sehr einfach durchführen. Listing 3 zeigt einige der erlaubten Möglichkeiten.
Listing 3: Polynom-Algebra
import numpy as <span class="hljs-built_in">np</span> <br/>from numpy.polynomial import Polynomial <br/><br/>p = Polynomial([<span class="hljs-number">6</span>, -<span class="hljs-number">5</span>, <span class="hljs-number">1</span>]) # <span class="hljs-number">6</span>-5x+x**<span class="hljs-number">2</span> <br/>q = Polynomial([<span class="hljs-number">2</span>, -<span class="hljs-number">3</span>]) # <span class="hljs-number">2</span>-3x <br/><br/># Additon, Subtraktion <span class="hljs-literal">und</span> Multiplikation <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Add:"</span>, p + q) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Sub:"</span>, p - q) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Mul:"</span>, p * q) <br/><br/># Multiplikation <span class="hljs-literal">und</span> Division mit Skalar <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"MulSkalar:"</span>, q * <span class="hljs-number">3.0</span>) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"DivSkalar:"</span>, q / <span class="hljs-number">2.0</span>) <br/><br/># Potenzierung <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Hoch 3:"</span>, p ** <span class="hljs-number">3</span>) <br/><br/># Division <br/><span class="hljs-built_in">quotient</span>, <span class="hljs-built_in">rest</span> = divmod(p, q) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Div: "</span>, <span class="hljs-built_in">quotient</span>) <br/><span class="hljs-built_in">print</span> (<span class="hljs-string">"Rest:"</span>, <span class="hljs-built_in">rest</span>)
Zunächst werden zwei Polynome p und q definiert. Die Algebra-Funktionen werden in Python zum großen Teil mithilfe von überladenen Standardoperatoren ermöglicht. So werden Addition, Subtraktion, Multiplikation und Potenzierung mit den Operationen +, -, * und ** ausgeführt. Außerdem gibt es die Möglichkeit, ein Polynom mit einem skalaren Wert zu verarbeiten.Interessanter ist jedoch die Division zweier Polynome. Das haben wir alle irgendwann mal in der Schule gelernt! Aber Hand aufs Herz: Wer weiß noch, wie es geht?Nun ist die Sache endgültig erledigt, denn Python übernimmt das für uns. Im letzten Teil von Listing 3 wird das Polynom p durch das Polynom q dividiert. Das Ergebnis ist ein neues Polynom und zusätzlich ein übrig bleibender Restwert. Führen Sie das Programm aus Listing 3 aus, so erhalten Sie die folgenden Ausgaben:
Add: poly([ 8. <span class="hljs-string">-8</span>. 1.])
Sub: poly([ 4. <span class="hljs-string">-2</span>. 1.])
Mul: poly([ 12. <span class="hljs-string">-28</span>. 17. <span class="hljs-string">-3</span>.])
MulSkalar: poly([ 6. <span class="hljs-string">-9</span>.])
DivSkalar: poly([ 1. <span class="hljs-string">-1</span>.5])
Hoch 3: poly([ 216. <span class="hljs-string">-540</span>. 558. <span class="hljs-string">-305</span>. 93.
<span class="hljs-string">-15</span>. 1.])
Div: poly([ 1.44444444 <span class="hljs-string">-0</span>.33333333])
Rest: poly([ 3.11111111])
Python gibt ein Polynom als poly(...)-Objekt aus, um es von der normalen Listenausgabe zu unterscheiden.Die Multiplikation der beiden Polynome p und q ergibt als Ergebnis ein neues poly-Objekt – das in der dritten Zeile in der Ausgabe oben:
p * q = <span class="hljs-number">12</span> – <span class="hljs-number">28</span>x + <span class="hljs-number">17</span>x<sup>2</sup> – <span class="hljs-number">3</span>x<sup>3</sup>
Wie wird nun der letzte Teil der Ausgabe, also das Ergebnis der Polynomdivision, interpretiert? Zunächst ist das Ergebnis von p / q wieder ein Polynom, und zwar 1.4444444 – 0.3333333x, oder etwas lesbarer: 13/9 – x/3. Der Rest der Polynomdivision ist 3.11111111 oder als Bruch: 28/9.
Die Nullstellen in Polynomen bestimmen
Und es wird noch spannender! Nun sollen die Nullstellen eines Polynoms berechnet werden. Dazu kommt das Produkt der beiden Polynome p * q aus Listing 3 zum Einsatz. Um einen Überblick über diese Funktion zu bekommen, soll sie zunächst einmal als Grafik ausgegeben werden (Listing 4).Listing 4: Polynom für die Nullstellensuche
<span class="hljs-keyword">import</span> numpy <span class="hljs-keyword">as</span> np <br/>from numpy.polynomial <span class="hljs-keyword">import</span> Polynomial <br/><span class="hljs-keyword">import</span> matplotlib.pylab <span class="hljs-keyword">as</span> plt <br/><br/>p = Polynomial([6, -5, 1]) # 6-5x+x**2 <br/>q = Polynomial([2, -3]) # 2-3x <br/><br/>x = np.linspace(0, 4, 101) <br/>y = (p*q)(x) <br/><br/>plt.plot(x,y) <br/>plt.grid() <br/>plt.show()
Hier wird eine recht interessante Syntax für die Berechnung der y-Werte benutzt: y = (p*q)(x). Es wird also zunächst die Polynommultiplikation p*q durchgeführt, dann das Ergebnis-Objekt mit dem Array der x-Werte gefüttert. Das liefert wiederum das Array der y-Werte, die dann als Grafik dargestellt werden, siehe Bild 1.
Wie Sie in der grafischen Darstellung erkennen können, hat das Polynom drei Nullstellen, was sich auch mit den Erwartungen deckt. Nun kann das Programm aus Listing 4 um zwei Zeilen erweitert werden, welche die Schnittstellen des Polynoms mit der x-Achse berechnen:
<span class="hljs-attr">nullstellen</span> = (p*q).roots()
print (<span class="hljs-string">"Nullst.:"</span>, nullstellen)
Als Ergebnis erhält man die Grafik und die drei Nullstellen: 0.6666667, 2.0 und 3.0.
Ableitungen und Integrale von Polynomen
Eigentlich sind Polynome leicht abzuleiten (differenzieren) oder zu integrieren. Dafür gibt es einige einfache Regeln, die angewendet werden müssen. Aber warum selbst rechnen, wenn es auch mit Python geht?Das Listing 4 (Nullstellen) soll nun noch um die Berechnung und die Darstellung der Ableitung und des Integrals erweitert werden. Listing 5 zeigt den kompletten Code dazu.Listing 5: Integral und Ableitung
<span class="hljs-built_in">import</span> numpy as np <br/>from numpy.polynomial <span class="hljs-built_in">import</span> Polynomial <br/><span class="hljs-built_in">import</span> matplotlib.pylab as plt <br/><br/><span class="hljs-attr">p</span> = Polynomial([<span class="hljs-number">6</span>, -<span class="hljs-number">5</span>, <span class="hljs-number">1</span>]) <br/><span class="hljs-attr">q</span> = Polynomial([<span class="hljs-number">2</span>, -<span class="hljs-number">3</span>]) <br/><span class="hljs-attr">x</span> = np.linspace(<span class="hljs-number">0</span>, <span class="hljs-number">4</span>, <span class="hljs-number">101</span>) <br/><span class="hljs-attr">y</span> = (p*q)(x) <br/><span class="hljs-attr">p1</span> = (p*q).deriv() <br/><span class="hljs-attr">p2</span> = (p*q).integ() <br/><br/>plt.plot(x,y) <br/>plt.plot(x,p1(x), <span class="hljs-attr">color=’k’)</span> <br/>plt.plot(x,p2(x), <span class="hljs-attr">color=’r’)</span> <br/>plt.grid() <br/>plt.show() <br/><br/><span class="hljs-attr">nullstellen</span> = (p*q).roots() <br/>print (<span class="hljs-string">"Nullst.:"</span>, nullstellen) <br/>print (<span class="hljs-string">"Ableitung:"</span>, p1) <br/>print (<span class="hljs-string">"Integral: "</span>, p2)
Die Berechnung der ersten Ableitung wird durch den Aufruf der Funktion deriv ermöglicht. Direkt danach wird die Integration durchgeführt. In beiden Fällen wird als Ergebnis wieder ein Polynomial-Objekt geliefert (hier: p1 und p2). Diese Polynome werden dann in die Grafik eingezeichnet (Funktion: blau, Ableitung: schwarz und Integral: rot). Nach den Nullstellen werden die Ergebnispolynome in der bekannten Schreibweise als poly-Objekte ausgegeben:
Ableitung: poly([<span class="hljs-string">-28</span>. 34. <span class="hljs-string">-9</span>.])
Integral: poly([ 0. 12. <span class="hljs-string">-14</span>. 5.66666667 <span class="hljs-string">-0</span>.75 ]) Lineare Algebra
In diesem Abschnitt geht es noch einmal um grundlegende Matrix-Operationen, die teilweise schon in [1] beschrieben wurden. Eine der wichtigsten Operationen mit Matrizen ist die Matrixmultiplikation, die in den Naturwissenschaften immer wieder benutzt wird. Aus diesem Grund sollte die Funktion für die Multiplikation zweier Matrizen auch möglichst schnell arbeiten.Eine Skalar-Multiplikation wird mit einem numpy-Array über den normalen Multiplikationsoperator ausgeführt:
import numpy as np
a = np.array(<span class="hljs-string">[[2.5, 1.0], [0.5, 3.0]]</span>)
<span class="hljs-built_in">print</span> (a)
b = a * <span class="hljs-number">3.0</span>; # Skalarmultiplikation
<span class="hljs-built_in">print</span> (b)↓
In diesem Fall wird einfach jedes Array-Element in a mit 3.0 multipliziert. Das Ergebnis-Array hat die gleiche Größe wie das Ausgangs-Array.Die Multiplikation zweier Matrizen wird mit np.dot oder np.matmul ausgeführt. Das Beispiel in Listing 6 zeigt drei mögliche Implementierungen: Zunächst eine sehr einfache Eigenimplementierung mit den drei üblichen, ineinander geschachtelten Schleifen. Danach kommen die einfachen Aufrufe der Funktionen im numpy-Modul.
Listing 6: Matrizen multiplizieren
import numpy as <span class="hljs-built_in">np</span> <br/><br/>a = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">array</span>([[<span class="hljs-number">2.5</span>, <span class="hljs-number">1.0</span>], [<span class="hljs-number">0.5</span>, <span class="hljs-number">3.0</span>]]) <br/>b = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">array</span>([[<span class="hljs-number">1.5</span>, <span class="hljs-number">2.0</span>], [<span class="hljs-number">1.5</span>, <span class="hljs-number">2.5</span>]]) <br/><br/># Eigener, sehr einfacher Algorithmus <br/>c = <span class="hljs-built_in">np</span>.zeros([<span class="hljs-number">2</span>, <span class="hljs-number">2</span>]) <br/><span class="hljs-keyword">for</span> i <span class="hljs-keyword">in</span> <span class="hljs-built_in">range</span>(<span class="hljs-number">2</span>): <br/> <span class="hljs-keyword">for</span> j <span class="hljs-keyword">in</span> <span class="hljs-built_in">range</span>(<span class="hljs-number">2</span>): <br/> <span class="hljs-keyword">for</span> k <span class="hljs-keyword">in</span> <span class="hljs-built_in">range</span>(<span class="hljs-number">2</span>): <br/> c[i, j] += a[i, k] * b[k, j] <br/><br/><span class="hljs-built_in">print</span> (c) <br/><br/># Eingebauter Algorithmus <br/>d = <span class="hljs-built_in">np</span>.dot(a, b) <br/><span class="hljs-built_in">print</span> (d) <br/><br/>e = <span class="hljs-built_in">np</span>.matmul(a, b) <br/><span class="hljs-built_in">print</span> (e)
Alle drei Implementierungen liefern wie erwartet das gleiche Ergebnis-Array:
[[ <span class="hljs-number">5.25</span> <span class="hljs-number">7.5</span> ]
[ <span class="hljs-number">5.25</span> <span class="hljs-number">8.5</span> ]]
Für zweidimensionale Arrays liefert die Funktion np.dot das gleiche Ergebnis wie die Funktion np.matmul. Bei eindimensionalen Vektoren wird das sogenannte innere Produkt ermittelt. Bei Arrays mit mehr als drei Dimensionen werden die Produkte über die beiden letzten Achsen summiert.
Da die Matrixmultiplikation für alle Bereiche der Naturwissenschaften und Technik enorm wichtig ist, sollen hier noch einmal einige Performance-Vergleiche durchgeführt werden. Grundsätzlich sollten Sie beim Betrachten der Ergebnisse allerdings mehrere Dinge im Hinterkopf haben:
- Man kann unterschiedlich leistungsfähige Algorithmen implementieren.
- Es werden unterschiedliche Sprachen und Laufzeitumgebungen verglichen.
- Man kann Parallelprogrammierung einsetzen.
Der Python-Code für das Beispiel befindet sich in der Datei CodeZuBild3.py im Download zu diesem Artikel und wird hier im Text nicht aufgelistet, da er Listing 6 stark ähnelt. Die erforderlichen großen Arrays wurden einfach mit Zufallszahlen initialisiert.Alle Spalten in Bild 3, die mit (1) gekennzeichnet sind, laufen in einer parallelen Umgebung mit vier Threads. Man kann erkennen, dass die normalen Rechenoperationen in Python relativ langsam sind, die entsprechenden Bibliotheksfunktionen jedoch eine sehr gute Performance zeigen. Die .NET-Sprache C# liegt etwa in der Mitte.Es stellt sich die Frage, warum die Python-Funktionen np.dot und np.matmul so schnell sind? Die Antwort ist einfach: Optimierung! Hier kommen viele Möglichkeiten zum Einsatz: Loop Unrolling, optimale Cache-Ausnutzung, richtige Indizierung der Arrays, Pallelisierung mit SSE und Threads oder spezielle Prozessorbefehle (wie zum Beispiel FusedMultiplyAndAdd). Es sollte klar sein, dass die Eigenimplementierung der Matrixmultiplikation (drei Schleifen) die einfachste und leider auch langsamste Variante ist. Dies gilt insbesondere für sehr große Arrays.
Gleichungssysteme lösen
Ein anderer wichtiger Bereich der linearen Algebra ist die Lösung von linearen Gleichungssystemen mit mehreren Unbekannten. Vielleicht erinnern Sie sich noch daran, wie Sie in der Schule Gleichungssysteme mit drei Unbekannten gelöst haben – mit viel Schweiß auf der Stirn!Das soll nun anders werden. Hierzu wird die Funktion np.linalg.solve eingesetzt. Ein Gleichungssystem mit drei Unbekannten sieht allgemein so aus:
<span class="hljs-name">m</span><sub><span class="hljs-name">11</span></sub> x<sub><span class="hljs-number">1</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">12</span></sub> x<sub><span class="hljs-number">2</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">13</span></sub> x<sub><span class="hljs-number">3</span></sub> = b<sub><span class="hljs-number">1</span></sub>
<span class="hljs-name">m</span><sub><span class="hljs-name">21</span></sub> x<sub><span class="hljs-number">1</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">22</span></sub> x<sub><span class="hljs-number">2</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">23</span></sub> x<sub><span class="hljs-number">2</span></sub> = b<sub><span class="hljs-number">2</span></sub>
<span class="hljs-name">m</span><sub><span class="hljs-name">31</span></sub> x<sub><span class="hljs-number">1</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">32</span></sub> x<sub><span class="hljs-number">2</span></sub> + <span class="hljs-name">m</span><sub><span class="hljs-name">33</span></sub> x<sub><span class="hljs-number">3</span></sub> = b<sub><span class="hljs-number">3</span></sub>
Diese drei Gleichungen können auch in einer Matrixschreibweise geschrieben werden: M x = b. Hierbei ist M die Matrix der Koeffizienten. Diese Matrix und der Vektor b müssen der Funktion np.linalg.solve zur Verfügung gestellt werden, um daraus den Lösungsvektor x zu ermitteln. Listing 7 zeigt ein Beispiel mit den drei folgenden Gleichungen:
Listing 7: Lineares Gleichungssystem lösen
import numpy as np <br/><br/># Arrays definieren <br/>M = np.array([[<span class="hljs-number">3</span>, <span class="hljs-number">-2</span>, <span class="hljs-number">0</span>], [<span class="hljs-number">-2</span>, <span class="hljs-number">1</span>, <span class="hljs-number">-3</span>], [<span class="hljs-number">4</span>, <span class="hljs-number">6</span>, <span class="hljs-number">1</span>]]) <br/>b = np.array([<span class="hljs-number">8</span>, <span class="hljs-number">-20</span>, <span class="hljs-number">7</span>]) <br/><br/># Gleichungssystem lösen <br/>x = np.linalg.solve(M, b) <br/><br/>print (x)
<span class="hljs-number">3</span>x – <span class="hljs-number">2</span>y = <span class="hljs-number">8</span>
<span class="hljs-number">-2</span>x + y - <span class="hljs-number">3</span>z = <span class="hljs-number">-20</span>
<span class="hljs-number">4</span>x + <span class="hljs-number">6</span>y + z = <span class="hljs-number">7</span>
Das wirklich einfache Python-Programm aus Listing 7 liefert die drei Unbekannten x, y und z als Lösungsvektor (Array):
[ <span class="hljs-number">2.</span> <span class="hljs-number">-1.</span> <span class="hljs-number">5.</span>] Interpolation
Als Nächstes soll das Modul scipy zum Einsatz kommen. Unter vielen anderen Möglichkeiten findet man hier auch die Definitionen verschiedener naturwissenschaftlicher Konstanten, und zwar mit dem Wert, der Maßeinheit und dem jeweiligen Messfehler.Im folgenden Beispiel sollen die Interpolationsmöglichkeiten (mit einer Variablen) im Modul scipy genauer betrachtet werden. Eine Interpolation von mehreren diskreten Datenwerten kann auf unterschiedliche Weisen ausgeführt werden und liefert dann auch voneinander abweichende Ergebnisse. Im gezeigten Beispiel (Listing 8) werden zwei Arrays angelegt, welche die gemessenen Werte (tmess und ymess) enthalten. Danach werden diese Daten als Grafik dargestellt. Nun kann die Interpolation (Funktion interp1d) ausgeführt werden. Es müssen die Daten-Arrays und die gewünschte Interpolationsart angegeben werden.Listing 8: Interpolation
<span class="hljs-keyword">import</span> numpy <span class="hljs-keyword">as</span> np <br/><span class="hljs-keyword">import</span> matplotlib.pylab <span class="hljs-keyword">as</span> plt <br/><span class="hljs-keyword">from</span> scipy.interpolate <span class="hljs-keyword">import</span> interp1d <br/><br/># Messwerte: <br/>tmess = np.array([<span class="hljs-number">0.0</span>, <span class="hljs-number">2.0</span>, <span class="hljs-number">4.0</span>, <span class="hljs-number">6.0</span>, <span class="hljs-number">8.0</span>, <span class="hljs-number">10.0</span>]) <br/>ymess = np.array([<span class="hljs-number">0.4</span>, <span class="hljs-number">25.0</span>, <span class="hljs-number">142.5</span>, <span class="hljs-number">155.6</span>, <span class="hljs-number">152.3</span>,<br/> <span class="hljs-number">112.2</span>]) <br/><br/># Darstellung als Punkte <br/>plt.plot(tmess, ymess, <span class="hljs-string">'o'</span>) <br/>plt.grid() <br/><br/># Interpolation <br/># kind=linear, quadratic, cubic, nearest... <br/>f_inter = interp1d(tmess, ymess, kind=<span class="hljs-string">'linear'</span>) <br/><br/># Berechnung und Darstellung der <br/># interpolierten Werte <br/>x = np.linspace(<span class="hljs-number">0.0</span>, <span class="hljs-number">10.0</span>, <span class="hljs-number">100</span>) <br/>plt.plot(x, f_inter(x)) <br/>plt.show() <br/><br/># Berechnung bestimmter interpolierter Werte <br/>print(<span class="hljs-string">"X = 3.0, Y ="</span>, f_inter(<span class="hljs-number">3.0</span>)) <br/>print(<span class="hljs-string">"X = 5.0, Y ="</span>, f_inter(<span class="hljs-number">5.0</span>))
benutzt werden. Das bedeutet, dass der Teil zwischen zwei Messpunkten jeweils durch eine Gerade interpoliert wird. Als Ergebnis erhält man eine Funktion, die eingesetzt werden kann, um beliebige Werte zwischen den diskreten Messwerten zu berechnen. Im Beispielcode wird die Ergebnisfunktion ebenfalls gezeichnet, zu sehen ist sie in Bild 4. Schließlich werden zwei interpolierte Werte an den Stellen 3.0 und 5.0 berechnet und ausgegeben.
Die Funktion interp1d aus dem Modul scipy ermittelt ein Rückgabeobjekt, das die interpolierte Funktion abbildet. Dieses Objekt kann also benutzt werden, um beliebige Interpolationswerte auszurechnen. Man muss jedoch beachten, dass die berechneten Werte im Definitionsbereich der interpolierten Daten liegen. Im gezeigten Beispiel müssen die x-Werte also zwischen 0.0 und 10.0 liegen.Die Bibliothek stellt unterschiedliche Interpolationsarten zur Verfügung. Die in Listing 8 angewendete Art linear verbindet die einzelnen diskreten Datenpunkte durch einfache Geradensegmente. Die Interpolationsart quadratic erzeugt quadratische, cubic erstellt kubische Spline-Funktionen, welche die Daten optimal abbilden.
Symbolische Mathematik mit sympy
Das Python-Modul sympy enthält viele Klassen, die symbolische Mathematik ermöglichen. Das bedeutet vereinfacht ausgedrückt, dass man direkt mit mathematischen Formeln rechnen kann. Sie können also einen mathematischen Ausdruck umstellen oder vereinfachen. Als Ergebnis erhalten Sie einen neuen mathematischen Ausdruck.Das hört sich ja sehr gut an! Wie oft haben wir in der Schule vor mathematischen Ausdrücken gesessen und mussten diese vereinfachen oder umstellen? Auch das wird nun dank Python und der Bibliothek sympy wesentlich einfacher. Hier ein einführendes Beispiel:
import sympy
# Ausdruck definieren
expr = <span class="hljs-number">2</span>*(<span class="hljs-keyword">x</span>**<span class="hljs-number">2</span> - <span class="hljs-keyword">x</span>) - <span class="hljs-keyword">x</span>*(<span class="hljs-keyword">x</span> + <span class="hljs-number">2</span>)
# Ausdruck vereinfachen
result = sympy.simplify(expr)
print (<span class="hljs-string">"Ausdruck: "</span>, expr)
print (<span class="hljs-string">"Vereinfacht:"</span>, result)
Zunächst wird das Modul sympy importiert. Dann definiert man einen mathematischen Ausdruck, der anschließend mit der Funktion simplify vereinfacht wird. Das Ergebnis lautet im vorliegenden Fall so:
x*(<span class="hljs-name">x</span> - <span class="hljs-number">4</span>)
Das sieht schon einmal sehr gut aus, aber leider ist es nicht in allen Fällen so einfach. Der folgende Code soll den Ausdruck sqrt(x**2) vereinfachen. Die Quadratwurzel einer quadrierten Zahl sollte die Zahl selbst ergeben. Als Ergebnis erwartet man also x. Der erste Teil des unten gezeigten Programmcodes verändert den Ausdruck aber gar nicht.Erst im zweiten Teil funktioniert alles wie gewünscht. Dort wurde das Symbol y so definiert, dass es immer einen positiven Wert enthalten soll. Lediglich unter dieser Bedingung ist die Wurzel-Funktion nämlich erlaubt und das Ergebnis ist dann x.
<span class="hljs-keyword">import</span> sympy
x = sympy.Symbol(<span class="hljs-string">"x"</span>)
<span class="hljs-built_in">print</span> (sympy.<span class="hljs-built_in">sqrt</span>(x**<span class="hljs-number">2</span>))
y = sympy.Symbol(<span class="hljs-string">"y"</span>, positive=<span class="hljs-literal">True</span>)
<span class="hljs-built_in">print</span> (sympy.<span class="hljs-built_in">sqrt</span>(y**<span class="hljs-number">2</span>))
Zum Abschluss möchte ich ein kleines (ganz einfaches) Mathematikprogramm erstellen, das mit Ausdrücken arbeiten kann. Es soll einen Ausdruck vereinfachen, ausmultiplizieren, faktorisieren, differenzieren oder integrieren können, siehe Listing 9. Die Benutzerschnittstelle wurde dabei bewusst einfach gehalten, da es hier nur um die Benutzung der Funktionalitäten im Modul sympy geht. Am Anfang des Programms wird das Symbol x explizit definiert. Dadurch vermeidet man im folgenden Code die Warnung, dass x undefiniert sei. Danach wird der mathematische Ausdruck als normale String-Variable eingelesen und die Rechenoperation abgefragt.
Listing 9: Einfaches Mathematik-Programm
<span class="hljs-keyword">import</span> sympy <br/><br/>x = sympy.<span class="hljs-type">Symbol</span>(<span class="hljs-string">"x"</span>) <br/><span class="hljs-built_in">expr</span> = input(<span class="hljs-string">"Ausdruck f(x): "</span>) <br/><br/>print (<span class="hljs-string">"Was soll gemacht werden:"</span>) <br/>print (<span class="hljs-string">" (v)ereinfachen"</span>) <br/>print (<span class="hljs-string">" (a)usmultiplizieren"</span>) <br/>print (<span class="hljs-string">" (f)aktorisieren"</span>) <br/>print (<span class="hljs-string">" nach x (d)ifferenzieren"</span>) <br/>print (<span class="hljs-string">" über x (i)ntegrieren"</span>) <br/>todo = input(<span class="hljs-string">"v, a, f, d oder i: "</span>) <br/><br/><span class="hljs-keyword">if</span>(todo == 'v'): <br/> <span class="hljs-literal">result</span> = sympy.simplify(<span class="hljs-built_in">expr</span>) <br/><span class="hljs-keyword">elif</span>(todo == 'a'): <br/> <span class="hljs-literal">result</span> = sympy.expand(<span class="hljs-built_in">expr</span>) <br/><span class="hljs-keyword">elif</span>(todo == 'f'): <br/> <span class="hljs-literal">result</span> = sympy.factor(<span class="hljs-built_in">expr</span>) <br/><span class="hljs-keyword">elif</span>(todo == 'd'): <br/> <span class="hljs-literal">result</span> = sympy.diff(<span class="hljs-built_in">expr</span>, x) <br/><span class="hljs-keyword">elif</span>(todo == 'i'): <br/> <span class="hljs-literal">result</span> = sympy.integrate(<span class="hljs-built_in">expr</span>, x) <br/> <br/>print (<span class="hljs-string">"Ergebnis:"</span>, <span class="hljs-literal">result</span>)
Mithilfe mehrerer if-Befehle wird nun die gewünschte Operation auf den Ausdruck angewendet und als Ergebnis ausgegeben. Ein Beispiellauf des Programms könnte wie folgt aussehen:
<span class="hljs-selector-tag">Ausdruck</span> <span class="hljs-selector-tag">f</span>(x): <span class="hljs-selector-tag">2</span>*<span class="hljs-selector-tag">exp</span>(-<span class="hljs-number">0.5</span>*x)/<span class="hljs-selector-tag">x</span>
<span class="hljs-selector-tag">Was</span> <span class="hljs-selector-tag">soll</span> <span class="hljs-selector-tag">gemacht</span> <span class="hljs-selector-tag">werden</span>:
(v)<span class="hljs-selector-tag">ereinfachen</span>
(a)<span class="hljs-selector-tag">usmultiplizieren</span>
(f)<span class="hljs-selector-tag">aktorisieren</span>
<span class="hljs-selector-tag">nach</span> <span class="hljs-selector-tag">x</span> (d)<span class="hljs-selector-tag">ifferenzieren</span>
ü<span class="hljs-selector-tag">ber</span> <span class="hljs-selector-tag">x</span> (i)<span class="hljs-selector-tag">ntegrieren</span>
<span class="hljs-selector-tag">v</span>, <span class="hljs-selector-tag">a</span>, <span class="hljs-selector-tag">f</span>, <span class="hljs-selector-tag">d</span> <span class="hljs-selector-tag">oder</span> <span class="hljs-selector-tag">i</span>: <span class="hljs-selector-tag">d</span>
<span class="hljs-selector-tag">Ergebnis</span>: <span class="hljs-selector-tag">-1</span><span class="hljs-selector-class">.0</span>*<span class="hljs-selector-tag">exp</span>(-<span class="hljs-number">0.5</span>*x)/<span class="hljs-selector-tag">x</span> <span class="hljs-selector-tag">-</span> <span class="hljs-selector-tag">2</span>*<span class="hljs-selector-tag">exp</span>(-<span class="hljs-number">0.5</span>*x)/<span class="hljs-selector-tag">x</span>**<span class="hljs-selector-tag">2</span> Zusammenfassung
Das ist erwartungsgemäß noch längst nicht alles, was Python in Sachen Mathematik draufhat. Die Möglichkeiten von numpy, scipy und sympy konnten in diesem zweiteiligen Artikel gerade einmal angerissen werden. Viele Klassen, die Verfahren für die höhere Mathematik bereitstellen, wurden hier nicht angesprochen, da sie nur in ganz speziellen Fällen benötigt werden. Einige Erkenntnisse sollen aber noch einmal herausgestellt werden:- Python ist bei mathematischen Berechnungen nicht langsam, sondern besonders flott.
- Python enthält sehr umfangreiche Mathe-Bibliotheken.
- Die Benutzung der Klassen ist sehr einfach.
- Viele Programmierfehler können durch die Benutzung der Bibliotheken vermieden werden.
- Letztendlich kann man eine Menge Zeit sparen.
Fussnoten
- Bernd Marquard, Rechnen mit Python, Teil 1, Arrays mit Schleife, dotnetpro 12/2017, Seite 64 ff., http://www.dotnetpro.de/A1712Rechnen
