18. Jun 2018
Lesedauer 13 Min.
Python kann noch mehr
Mathematik mit Python, Teil 3
Wir loten die Möglichkeiten der numerischen Klassen in der Python-Bibliothek tiefer aus: Erläuterungen und Beispiele zu den Themen Optimierung, Fitting, Integration und Graphen.
Bild: Shutterstock / Alexander Supertramp
Auch in diesem dritten Artikel zu den mathematischen Möglichkeiten von Python wird es wieder um Performance gehen. Darüber hinaus soll dieses Mal aber auch ein Blick auf die komplexeren mathematischen Funktionalitäten der Python-Bibliotheken numpy, scipy, sympy und matplotlib geworfen werden. Aber keine Angst: Es bleibt alles noch auf einem verständlichen Niveau.Wenn Sie die Beispiele aus diesem Artikel und den beiden vorangegangenen Beiträgen [1][2] ausprobieren wollen, sollten Sie unbedingt die neuesten Versionen der Bibliotheken numpy, scipy, sympy und matplotlib installieren.
Zufallszahlen
Wenn man es ganz genau nimmt, ist die Berechnung von „guten“ Zufallszahlen gar nicht so einfach, denn zufällige Zahlen lassen sich nun mal nicht „berechnen“! Man erhält normalerweise sogenannte Pseudozufallszahlen, mit denen man dann die gewünschten Berechnungen ausführen kann.Um Zufallszahlen zu berechnen, können Sie die Funktionen random.rand() und random.randint(von, bis) nutzen. Diese Funktionen liefern Fließkommazahlen zwischen 0,0 und kleiner 1,0 oder Ganzzahlen im angegebenen Bereich.Hier ist allerdings noch eine wichtige Sache zu berücksichtigen. Der Zufallszahlengenerator wird zufällig initialisiert, das heißt, mehrere Programmläufe erzeugen unterschiedliche Zahlensequenzen:<span class="hljs-keyword">import</span> numpy <span class="hljs-keyword">as</span> np
# Immer unterschiedliche Zahlen
z1 = np.random.randint(0, 1000)
print (z1)
z2 = np.random.randint(0, 1000)
print (z2)
Für Testzwecke ist es aber in vielen Fällen sinnvoll, wiederholbare Zahlenfolgen zu generieren. Um dies zu erreichen, können Sie den Zufallszahlengenerator über den Aufruf random.seed(start) initialisieren. Solange Sie den Parameter start nicht ändern, wird beim erneuten Programmstart immer wieder die gleiche Zahlenfolge erstellt:
import numpy as <span class="hljs-built_in">np</span>
# Immer die gleiche Zahlenfolge
<span class="hljs-built_in">np</span>.<span class="hljs-built_in">random</span>.seed(<span class="hljs-number">111</span>)
z1 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">random</span>.randint(<span class="hljs-number">0</span>, <span class="hljs-number">1000</span>)
<span class="hljs-built_in">print</span> (z1)
z2 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">random</span>.randint(<span class="hljs-number">0</span>, <span class="hljs-number">1000</span>)
<span class="hljs-built_in">print</span> (z2)
In Listing 1 wird die zufällige Verteilung der Würfe mit zwei Würfeln berechnet und grafisch dargestellt. Es wird zunächst ein Integer-Array angelegt, das 13 Elemente enthält, also von 0 bis einschließlich 12. Danach werden zweimal 10 000 Zufallszahlen erzeugt. Die Summe von jeweils zwei Zahlen wird nun benutzt, um den Zähler im jeweiligen Array-Element um eins zu erhöhen.
Listing 1: Zwei Würfel
import numpy as <span class="hljs-built_in">np</span> <br/>import matplotlib.pyplot as plt <br/><br/># Array für die Augenzahl von <span class="hljs-number">2</span> ... <span class="hljs-number">12</span> <br/>y = <span class="hljs-built_in">np</span>.zeros(<span class="hljs-number">13</span>, dtype=int) <br/><br/># Zwei Zufallszahlen berechnen <br/># Summe im Array aufaddieren <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">1</span>, <span class="hljs-number">10000</span>): <br/> w1 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">random</span>.randint(<span class="hljs-number">1</span>, <span class="hljs-number">7</span>, dtype=int) <br/> w2 = <span class="hljs-built_in">np</span>.<span class="hljs-built_in">random</span>.randint(<span class="hljs-number">1</span>, <span class="hljs-number">7</span>, dtype=int) <br/> y[w1+w2] += <span class="hljs-number">1</span> <br/><br/><span class="hljs-built_in">print</span> (y[<span class="hljs-number">2</span>:<span class="hljs-number">13</span>]) <br/><br/># Darstellung des Ergebnisses <br/>x = <span class="hljs-built_in">np</span>.linspace(<span class="hljs-number">0</span>, <span class="hljs-number">12</span>, <span class="hljs-number">13</span>, dtype=int) <br/>plt.bar(x[<span class="hljs-number">2</span>:<span class="hljs-number">13</span>], y[<span class="hljs-number">2</span>:<span class="hljs-number">13</span>]) <br/>plt.<span class="hljs-built_in">show</span>()
Nach der numerischen Ergebnisausgabe erzeugt der Aufruf von plt.bar(x, y) ein schlichtes Balkendiagramm. Das Integer-Array x enthält die x-Positionen der Balken und wird einfach mit der linspace-Funktion erstellt.Die Wahrscheinlichkeiten für die unterschiedlichen Augenkombinationen steigen linear an. Nach dem Maximum bei der Augenzahl sieben fällt die Wahrscheinlichkeit wieder linear ab. Für die Augenzahl zwei beträgt die Wahrscheinlichkeit 1/36, für drei 2/36, für vier 3/36 und so weiter.
Least Square Fits
Mit der „Methode der kleinsten Quadrate“ [3] lassen sich gegebene Datenpaare x und y (zum Beispiel Messwerte) an eine vorgegebene Funktion y = f(x) anpassen.Im Beispiel aus Tabelle 1 sind elf Messwerte-Paare vorgegeben.Tabelle 1: Messwerte
|
Eine grafische Darstellung dieser Messdaten kann sehr schnell mit den Plot-Funktionen aus der Bibliothek matplotlib erstellt werden (Bild 1).
Wie aus dieser Darstellung unschwer zu erkennen ist, haben die gemessenen Daten annähernd die Funktion einer Gauß-Kurve [4]. Solche Funktionen haben die allgemeine Form
f(x) = y0 + a * e-b * x * x
Es handelt sich also um eine Exponentialfunktion mit x2 im Exponenten und den drei unbekannten Konstanten y0, a und b. Genau diese Konstanten müssen bestimmt werden, um die vollständige Ausgleichsfunktion zu erhalten.Im ausgeglichenen Bereich kann man dann neue Werte problemlos interpolieren. Eine Extrapolation außerhalb des Bereichs sollte allerdings mit größter Vorsicht durchgeführt werden.Für die Form dieser Ausgleichsfunktion muss natürlich immer ein bisschen nachgedacht werden. Im einfachsten Fall möchte man vielleicht eine Gerade als Funktion benutzen:
# Lineare Funktion
def f(x, a, b):
return a + b * x
Natürlich sind auch quadratische Funktionen möglich:
# Quadratische Funktion
def f(x, a, b):
return a + b * x**2
In den meisten Fällen braucht man jedoch etwas kompliziertere Ausgleichsfunktionen. Hier ist also immer ein wenig „mathematische Kreativität“ angesagt.Bei einem „Least Square Fit“ werden die drei Konstanten der oben gezeigten Gauss-Funktion y0, a und b nun so bestimmt, dass die Summe der Quadrate aller Messwerte minus dem dazugehörigen berechneten Wert so klein wie möglich wird. Das hört sich zunächst einmal kompliziert an, lässt sich aber mit Python und der Bibliothek scipy ziemlich leicht programmieren (Listing 2).
Listing 2: Least Square Fit
import matplotlib.pyplot as plt <br/>from scipy import optimize <br/>import numpy as np <br/><br/># Die Funktion, welche die Daten abbilden soll <br/>def f(x, y0, a, b): <br/> return y0 + a * np.exp(-b * x**2) <br/><br/># Hilfsfunktionen: Differenz zwischen Soll und Ist <br/>def delta(coeff): <br/> return ydata - f(xdata, *coeff) <br/><br/># Messdaten <br/>xdata = np.linspace(0, 5, 11) <br/>ydata = [1.0, 0.905, 0.71, 0.552, 0.351, 0.303, 0.252, 0.23, 0.241, 0.22, 0.22] <br/><br/># Least Square Fit <br/>coeff_start=(1, 1, 1) <br/>coeff_opt, coeff_cov = optimize.leastsq(delta, coeff_start) <br/>print ("y0, a, b = ", coeff_opt) <br/><br/># Grafische Ausgabe <br/>fig, ax = plt.subplots() <br/>ax.scatter(xdata, ydata, color='k') <br/>ax.set_xlabel(r"$x$", fontsize=18) <br/>ax.set_ylabel(r"$f(x)$", fontsize=18) <br/>ax.grid() <br/>xxdata = np.linspace(0, 5, 101) <br/>ax.plot(xxdata, f(xxdata, *coeff_opt), 'r', linewidth=2) <br/>plt.savefig("Bild2.png")
Nach den diversen import-Befehlen müssen zwei Funktionen definiert werden. Im ersten Schritt wird die Funktion f(x) implementiert. Da die Daten in einer Gauß-Funktion approximiert werden sollen, wird hier die weiter oben gezeigte Exponentialfunktion mit den drei Konstanten y0, a und b programmiert.Die zweite Hilfsfunktion, die bereitgestellt werden muss, ist immer die gleiche. Hier wird einfach die Differenz zwischen dem Original-Datenwert und dem jeweils berechneten Wert aus der Funktion f(x, y0, a, b) ermittelt.Im Hauptteil des Programms werden zunächst die Messdaten in zwei Arrays xdata und ydata bereitgestellt. Danach wird das Tupel coeff_start deklariert, das die geschätzten Startwerte der Koeffizienten für die Approximation enthält. Mit diesem Koeffizientensatz beginnt die Berechnung. Die Funktion optimize.leastsq(...) aus der Bibliothek scipy verbessert anschließend diese Werte, bis das Quadrat der Unterschiede zwischen den Messwerten und berechneten Werten minimal wird.Die optimierten Koeffizienten werden dann als Zahlen (hier gerundet) ausgegeben:
y0 = 0.229 a = 0.758 b = 0.415
Die Konstante y0 gibt den Abstand der optimierten Kurve zur x-Achse an. Dabei ist a ein Maß für die Höhe und b die Breite der Gauß-Kurve. Eine größere Konstante b ergibt eine schmalere Kurve.Im letzten Teil des Programms werden die Messwerte und die berechneten Funktionswerte in einer Grafik dargestellt. Dabei wird einmal die Funktion scatter benutzt, um die diskreten Messdaten auszugeben. Die berechnete Funktion kann mit der plot-Funktion als geschlossene Linie dargestellt werden (Bild 2).
Optimierung
Können Sie sich noch an die erste Minmax-Aufgabe erinnern, die Sie im Mathematikunterricht gelöst haben? Ja, genau: das legendäre Dosenproblem! Die Schüler sollten eine Dose mit einem Liter Inhalt berechnen, die aber möglichst wenig Blech bei der Herstellung verbraucht. Dafür benötigte man die Differentialrechnung, also die erste Ableitung der Funktion für die Fläche.Nun können Sie die Lösung mit Python auf zwei Wegen erzielen. Zum einen können Sie analytisch an die Sache herangehen. Hierzu benutzen Sie die Funktionalität der Bibliothek sympy. Dieser Weg entspricht der Vorgehensweise in der Schule. Anderseits können Sie die Optimierungsfunktionen aus der Bibliothek scipy anwenden. Dadurch sparen Sie sich das Differenzieren, denn man benötigt nur die Funktion der Fläche in Abhängigkeit vom Radius r.Nach den import-Deklarationen wird in Listing 3 zunächst die Funktion programmiert, welche die Fläche nur in Abhängigkeit vom Radius r darstellt. Diese Funktion benötigen wir jedoch nicht bei der analytischen Lösung, die zuerst gezeigt werden soll. Die Formel in der Funktion wird folgendermaßen hergeleitet:Listing 3: Lösung des Dosenproblems
from scipy import optimize <br/>import numpy as np <br/>import sympy <br/><br/># Funktion flaeche(r) <br/>def f(r): <br/> return 2 * np.pi * r**2 + 2 / r <br/><br/># Analytische Lösung mit sympy <br/>r, h = sympy.symbols("r, h") <br/>flaeche = 2 * sympy.pi * r**2 + 2 * sympy.pi * r * h <br/>volumen = sympy.pi * r**2 * h <br/>h_r = sympy.solve(volumen - 1.0)[0] <br/>flaeche_r = flaeche.subs(h_r) <br/>rsol = sympy.solve(flaeche_r.diff(r))[0] <br/>print ("Analytische Methode:") <br/>print (rsol.evalf()) <br/>print (flaeche_r.subs(r, rsol).evalf()) <br/><br/># Numerische Lösung <br/>r_min = optimize.brent(f, brack=(0.1, 4)) <br/>print ("Numerische Methode:") <br/>print (r_min) <br/>print (f(r_min)) f(r,h) = 2 * p * r2 + 2 * p * r * h
v(r,h) = p * r2 * h
Da das Volumen v der Dose einen Liter betragen soll, kann man die zweite Gleichung umformen:
h(r) = 1.0 / (p * r2)
Nun ist in der letzten Gleichung die Höhe h nur vom Radius abhängig. Man kann dieses Ergebnis in die erste Gleichung einsetzen und erhält dann die Formel für die Fläche f, nur vom Radius r abhängig. Nach ein bisschen Kürzen ergibt sich:
f(r) = 2 * p * r2 + 2 / r
Diese Formel wird in der Funktion f(r) am Anfang des Python-Programms definiert.Zunächst möchte ich nun die analytische Lösung vorstellen. Da die Bibliothek sympy benutzt werden soll, müssen wir als Erstes die verwendeten Variablen r und h definieren. Jetzt werden die beiden wichtigen Formeln für die Oberfläche und das Volumen der Dose in symbolischer Form programmiert. Im nächsten Schritt wird die Volumen-Formel nach der Höhe h aufgelöst und in die Formel für die Oberfläche eingesetzt. In der Variablen flaeche_r befindet sich die letzte Gleichung aus der obigen Herleitung in symbolischer Form. Zur Überprüfung kann man die Formel mit print(flaeche_r) ausgeben.Nun muss die erste Ableitung ermittelt werden. Dies wird durch die Anwendung der diff-Funktion auf die Variable flaeche_r erreicht. Da das Minimum in der Flächenfunktion gesucht wird, muss die erste Ableitung an dieser Stelle null sein. Aus diesem Grund wird die solve-Funktion für die Ableitung aufgerufen. Man erhält also die Nullstelle der abgeleiteten Flächenfunktion.Dieses Ergebnis gibt den Radius (etwa 0,5419 Dezimeter, weil wir 1 Liter als Volumen eingesetzt haben, also 5,419 Zentimeter) der gesuchten optimalen Dose an. Die dazugehörige Höhe kann man mit der dritten Formel berechnen, und dies ergibt 10,838 Zentimeter.Nun kommt die numerische Lösung des Problems dran, die wesentlich einfacher ist. Es stehen verschiedene Optimierungsroutinen zur Verfügung. Hier soll die Funktion optimize.brent(...) benutzt werden. Die Parameter im Funktionsaufruf sind schnell erklärt: Zuerst wird die oben angelegte Flächenfunktion f(r) angegeben. Als zweiter Parameter wird der zu durchsuchende Bereich vorgegeben (hier 0.1 bis 4.0). Irgendwo in diesem Bereich sollte die Lösung liegen.Mit der angewendeten Lösungsvariante erhalten wir mit vier Stellen Genauigkeit ein identisches Ergebnis. Erst nach der achten Stelle ergeben sich dann jedoch geringe Unterschiede.Im gezeigten Beispiel heißt es jedoch etwas vorsichtig sein. Der Suchbereich wurde nämlich ganz bewusst mit dem Wert 0.1 begonnen, da der Wert 0.0 einen ZeroDivisionError ergeben würde.
Integration
Nur noch einmal zur Erinnerung: Unter Integration versteht man (im numerischen Sinn) die Bestimmung der Fläche unterhalb einer Kurve innerhalb von bestimmten Grenzen (vergleiche Bild 3).
Wenn Sie mit den Python-Bibliotheken sympy und scipy arbeiten, gibt es wiederum zwei Möglichkeiten, um Integrale zu lösen. Die erste Variante arbeitet analytisch und nutzt die Bibliothek sympy, die zweite Variante mit der scipy-Bibliothek löst das Problem numerisch. Die zweite Vorgehensweise bietet hierbei einen entscheidenden Vorteil: Man kann (fast) jede Funktion numerisch innerhalb vorgegebener Grenzen integrieren. Die Fläche unter der Kurve wird hier in sehr schmale Rechtecke unterteilt, deren Teilflächen sich einfach berechnen lassen und die dann zum Endergebnis aufaddiert werden.Eine analytische Lösung ist bei komplizierten Funktionen oft nur sehr schwer zu ermitteln, und in vielen Fällen geht es auch gar nicht.In den folgenden Beispielen werden darum beide Vorgehensweisen und ihre Grenzen vorgestellt.Zunächst soll das Integral der Funktion f(x) = –x2 + 5 in den Grenzen von x0 = –1 und x1 = 1.5 berechnet werden. In Bild 3 sind die entsprechenden Grenzlinien in die Funktion eingezeichnet.Am Anfang von Listing 4 wird zunächst die Funktion f(x) definiert. Diese Funktion wird vom numerischen Verfahren benutzt. Im ersten Teil des Programms wird die numerische Integration mit der integrate.quad-Methode aus der Bibliothek scipy ausgeführt. Als Rückgabewerte erhält man einmal die Fläche unter der Kurve von x = –1 bis x = 1.5 und andererseits eine Abschätzung des Rechenfehlers.
Listing 4: Integration der Funktion f(x) = -x2+5
from scipy import integrate <br/>import sympy <br/><br/>def f(x): <br/> return -x**2 + 5 <br/><br/># Numerische Integration <br/>val, err = integrate.quad(f, -1, 1.5) <br/>print ("Fläche (numerisch): ", val, " Fehler: ", err) <br/>print () <br/><br/># Analytische Integration <br/>X = sympy.symbols("x") <br/>expr = "-x**2 + 5" <br/>res = sympy.integrate(expr, X) <br/>print ("Integrierte Funktion: ", res) <br/><br/># Numerischen Wert in den Grenzen bestimmen <br/>start = res.subs(X, -1).evalf(15) <br/>ende = res.subs(X, 1.5).evalf(15) <br/>print ("Fläche (analytisch): ", ende - start)
Im zweiten Teil des Beispiels aus Listing 4 wird das analytische Verfahren benutzt. Zunächst wird das Symbol "x" deklariert. Danach wird in der Variablen expr der Ausdruck definiert, der integriert werden soll. Nun wird dieser Ausdruck in expr analytisch mit der sympy.integrate-Methode integriert. Das Ergebnis in der Variablen res ist hier wieder ein Ausdruck, der ausgegeben werden kann.Jetzt muss noch die Fläche unter der Kurve in den vorgegebenen Grenzen berechnet werden. Dies wird durch Einsetzen der jeweiligen x-Werte (Variablen ende und start im Listing) und anschließende Subtraktion erreicht.Beachten Sie bitte an dieser Stelle, dass die bei der Berechnung ermittelten Flächen auch negativ sein können. Dies hängt immer von den jeweils benutzten Grenzwerten in der x-Richtung ab.Im nächsten Beispiel soll die Kreiszahl Pi mithilfe der Integrationsverfahren berechnet werden. Dazu geht man von der Pythagoras-Formel r2 = x2 + y2 aus. Für den Radius r wird 1 für den Radius des Einheitskreises eingesetzt. Danach wird nach y aufgelöst. Das Ergebnis ist:
y = (1 – x2)½
Um die Zahl Pi zu berechnen, muss man nun einfach das Integral dieser Funktion in den Grenzen von x0 = 0 bis x1 = 1 berechnen und dieses Ergebnis mit dem Wert 4 multiplizieren, da ja nur die Fläche eines Viertelkreises berechnet wurde. Dies ist wiederum analytisch und numerisch möglich (Listing 5). Allerdings sind in diesem Beispiel die Integrationsgrenzen sehr genau zu beachten. Diese Grenzen dürfen auf keinen Fall kleiner als –1 und größer als +1 werden, da dann die Zahl in der Wurzel negativ wird und die Berechnung dann nicht mehr möglich ist.
Listing 5: Berechnung der Zahl Pi durch Integration
import numpy as np <br/>from scipy import integrate <br/>import sympy <br/><br/>def fpi(x): <br/> return np.sqrt(1 - x**2) <br/><br/># Numerische Integration <br/>valpi, err = integrate.quad(fpi, 0, 1) <br/>print ("Pi (numerisch): ", valpi * 4) <br/><br/># Analytische Integration <br/>X = sympy.symbols("x") <br/>expr = sympy.sqrt(1 - X**2) <br/>respi = sympy.integrate(expr, X) <br/>print ("Integrierte Funktion: ", respi) <br/>start = respi.subs(X, 0).evalf(16) <br/>ende = respi.subs(X, 1).evalf(16) <br/>print ("Pi (analytisch): ", 4 * (ende - start)) <br/><br/># Pi als Konstante aus numpy <br/>print("Pi (Konstante): ", np.pi)
Die Vorgehensweise ist im analytischen und im numerischen Fall identisch zum vorangegangenen Beispiel. Wenn man das Programm ausführt, erkennt man, dass die integrierte Funktion selbst bei der benutzten einfachen Formel schon ziemlich kompliziert ist. Das bedeutet, dass bei komplexen Funktionen eine Integration in der Regel mit dem numerischen Verfahren durchgeführt werden muss.Ein Vergleich der Ergebnisse liefert eine fast hundertprozentige Übereinstimmung zur vorgegebenen Konstante Pi aus der Bibliothek numpy.Beide Verfahren haben jedoch ihre Grenzen. Bei der Anwendung des analytischen Verfahrens muss das Integral gelöst werden können. Der sich daran anschließende Rechenaufwand ist ziemlich gering. Beim numerischen Verfahren kann der Rechenaufwand bei komplizierten Funktionen sehr groß werden. Außerdem muss sichergestellt werden, dass die Funktion zwischen den gegebenen Grenzwerten auch berechenbar ist.Bei der numerischen Integration kann man als Grenzwerte sogar -8 (-np.inf) und +8 (np.inf) benutzen.
def f(x):
return np.exp(-x**2)
val, err = integrate.quad(f, -np.inf, np.inf)
In solchen Fällen muss die Funktion jedoch im Unendlichen gegen null gehen, um eine numerische Integration zu ermöglichen. Das analytische Verfahren ist dann oft wesentlich leistungsfähiger.
Graphentheorie
Haben Sie sich auch schon einmal gefragt, welches der kürzeste Weg von „A“ nach „B“ ist? In den letzten Jahren vermutlich eher nicht, da heutzutage in jedem Smartphone ein Navigationssystem vorhanden ist. Aber wie funktioniert das eigentlich?Hier kommt die Python-Bibliothek networkx mit der Klasse Graph() zum Einsatz. Die Graphentheorie [5] ist eine sehr spannende Angelegenheit und die Vorgehensweise soll zunächst an einem ganz einfachen Beispiel gezeigt werden:import networkx as nx
g = nx.Graph()
g.add_nodes_from([1, 2, 3, 4, 5])
g.add_edges_from([(1, 2), (2, 3), (3, 4), (1, 5),
(4, 5)])
print (nx.shortest_path(g, 1, 4))
In dem kleinen Programm werden nach der Deklaration des Graph()-Objekts zuerst fünf Knotenpunkte mit der Funktion add_nodes_from definiert. Die Funktion nimmt als Parameter ein Array mit beliebigen Elementen an. Hier werden einfach Ziffern benutzt. Man kann aber beispielsweise auch Texte oder andere Datentypen verwenden. Danach müssen die Verbindungen zwischen den Knoten angegeben werden. Bild 4 zeigt die gewünschte Anordnung der Knoten und Verbindungen für dieses Beispiel.
Die Verbindungen werden anschließend mit der Funktion add_edges_from vorgegeben. Hier wird ein Array von Tupeln als Parameter übergeben.Schließlich soll der kürzeste Weg zwischen den Knoten (1) und (4) ermittelt werden. Hierbei ist zu beachten, dass in diesem Beispiel alle Kanten (oder Verbindungen) die gleiche Wertigkeit haben.Wenn Sie nun nach der kürzesten Verbindung suchen, bekommen Sie immer den Weg, der die minimale Anzahl von Knoten durchläuft. Die Funktion shortest_path(g, 1, 4) gibt als Ergebnis aus:
[1, 5, 4]
Anfangs- und Endknoten werden dabei im Ergebnis mit ausgegeben.Fürs Erste ist das doch schon ganz gut. Wenn Sie jedoch den Weg in einem realen Straßennetz suchen, ist die Angelegenheit ein klein wenig komplizierter, denn man benötigt für die einzelnen Verbindungen noch eine Gewichtung.Im einfachsten Fall kann man hierfür den Kilometerabstand zwischen den Kontenpunkten wählen. Diese Gewichtungsfaktoren können aber auch noch wesentlich erweitert werden: Straßenzustand, Höhenunterschied, Staugefahr, aktuelle Staulage, Straßentyp (Bundesstraße, Autobahn) oder Schwierigkeitgrad.In unserem einfachen Beispiel soll nur der Kilometerabstand zwischen den Knoten benutzt werden. Die entsprechenden Daten wurden in der Datei data.txt abgelegt. Die ersten Zeilen dieser Datei sehen folgendermaßen aus:
KI,HH,98 #Stadt1, Stadt2, km
HH,HL,70
HL,RO,122
RO,B,235
HH,B,289
HH,H,159
H,B,286
B,DD,193
...
Es handelt sich um eine kommaseparierte Datei, die sich mit Python sehr leicht einlesen lässt. Die Datei enthält übrigens keine Knotendefinitionen. Die Graph()-Klasse legt die erforderlichen Knoten automatisch an, wenn eine Verbindung definiert wird. Für die einzelnen Städte werden die Autokennzeichen benutzt. Der letzte Parameter in jeder Zeile ist der Abstand zwischen den Städten in Kilometern.In Listing 6 werden zunächst die Bibliotheken networkx und csv importiert. Danach wird das Graph()-Objekt initialisiert. Nun kann die Datei mit den Orten und den Entfernungen geöffnet werden. In diesem Fall liest man die Daten am besten mit der csv.reader-Funktion zeilenweise ein. Die drei Elemente einer Zeile stehen dann im Array row zur Verfügung und können direkt als Tupel mit der Funktion add_weighted_edges_from in den Graphen eingefügt werden. Nach dem Schließen der Datei werden zunächst alle Knoten des Graphen ausgegeben.
Listing 6: Die Suche des optimalen Weges
import networkx as nx <br/>import csv <br/><br/>g = nx.Graph() <br/><br/>infile = open('data.txt', 'r') <br/>for row in csv.reader(infile): <br/> a = (row[0], row[1], float(row[2])) <br/> g.add_weighted_edges_from([a]) <br/> <br/>infile.close() <br/><br/>print (g.nodes()) # Alle Nodes listen (Test) <br/><br/>start = input("Start-Ort: ") <br/>ziel = input("Ziel-Ort: ") <br/><br/>try: <br/> print ("---------------------------------------------") <br/><br/> p = nx.dijkstra_path(g, start.upper(), ziel.upper()) <br/> sum = nx.dijkstra_path_length(g, start.upper(), ziel.upper()) <br/><br/> print ("Route: ", p) <br/> print ("Distanz: ", sum, " Kilometer") <br/>except: <br/> print ("Fehler bei der Eingabe.") <br/><br/>print ("---------------------------------------------")
Nun werden Start- und Zielort als Autokennzeichen abgefragt. Da im Beispiel gewichtete Verbindungen zwischen den Knoten definiert wurden, wird nun der nach seinem Erfinder Edsger W. Dijkstra benannte Dijkstra-Algorithmus angewendet. Dieser Suchalgorithmus benutzt die vorgegebenen Gewichtungen, um einen Weg zu suchen, bei dem die Summe der Gewichtungswerte einen möglichst kleinen Wert ergibt. Die Suche wird im Beispiel in einem try-except-Block ausgeführt, falls die eingegebenen Städtekennzeichen nicht in der Liste vorhanden sind.Probieren Sie es einfach mal aus: Von Köln ('K') nach Berlin ('B') gibt das Programm aus:
--------------------------------
Route: ['K', 'DO', 'H', 'B']
Distanz: 598.0 Kilometer
--------------------------------
Das Ergebnis besteht aus zwei Teilen. Zunächst wird mit der Funktion nx.dijkstra_path(...) der Weg zum Ziel ermittelt. Die Parameter bei diesem Funktionsaufruf sind einfach nur das Graphen-Objekt sowie Start- und Zielort. Die eingegebenen Städtekürzel werden sicherheitshalber in Großbuchstaben umgewandelt.In unserem Beispiel erhalten wir als die Reiseroute: Köln, Dortmund, Hannover, Berlin. Im zweiten Aufruf (nx.dijkstra_path_length) wird die Gesamtlänge des optimalen Weges ermittelt. Die Angabe von 598 km passt ziemlich gut zum realen Straßennetz Deutschlands.
Fazit
Was schon im vorangegangenen Artikel gezeigt wurde, hat sich wieder bestätigt: Die vielen Python-Bibliotheken sind extrem leistungsfähig. Viele komplexe Berechnungen lassen sich oft mit wenigen Zeilen Programmcode umsetzen.Vergessen Sie nicht: Je weniger man programmieren muss, desto weniger Fehler macht man.Fussnoten
- Bernd Marquardt, Mathematik mit Python, Teil 1, Arrays mit Schleife, dotnetpro, 12/2017, Seite 64 ff., http://www.dotnetpro.de/A1712Rechnen
- Bernd Marquardt, Mathematik mit Python, Teil 2, Der Mathe-Turbo, dotnetpro 1/2018, Seite 78 ff., http://www.dotnetpro.de/A1801Rechnen
- Wikipedia: Methode der kleinsten Quadrate, https://de.wikipedia.org/wiki/Methode_der_kleinsten_Quadrate
- Wikipedia: Normalverteilung, https://de.wikipedia.org/wiki/Normalverteilung
- Wikipedia: Graphentheorie, https://de.wikipedia.org/wiki/Graphentheorie
