3. Lineare Algebra mit NumPy#
3.1. NumPy installieren#
In vielen mathematischen Anwendungen muss mit Vektoren und Matrizen gerechnet werden, beispielsweise bei der numerischen Berechnung von Integralen, Differentialgleichungen oder Problemen aus der Graphentheorie. Wir wollen uns in diesem Kapitel mit einer Python-Bibliothek beschäftigen, welche entsprechende Datentypen für Vektoren und Matrizen, sowie die üblichen Rechenoperationen für diese bereitstellt. Die Rede ist vom Paket NumPy.
Zunächst muss die NumPy-Bibliothek installiert werden. Conda erledigt dies mit dem Konsolenbefehl:
conda install numpy
Nun muss die Bibliothek in unser Python-Skript eingebunden werden. Wir könnten dies wie in Abschnitt Variablen und Datentypen mit der Zeile from numpy import * machen, was aber nicht empfehlenswert ist, da die NumPy-Bibliothek die Funktionen sin, cos, sqrt, etc. bereit stellt, welche die aus der math-Bibliothek überschreiben würden. Daher nutzen wir folgende Variante:
import numpy as np
Der Zusatz as np gibt nur an, dass wir die Bibliothek in Zukunft unter dem kürzeren Namen np und nicht unter dem langen numpy ansprechen können.
3.2. Arbeiten mit Vektoren#
3.2.1. Vektoren erzeugen#
Ein Vektor bzw. NumPy-Array kann mit der Funktion np.array(...) von einem beliebigen iterierbaren Objekt (Liste, Tupel, …), deren Elemente vom gleichen Typ sind, initialisiert werden:
a = np.array([1.,2.,3.])
b = np.array((6.,5.,4.))
print("a ist ein", type(a), "mit Wert", a)
print("b ist ein", type(b), "mit Wert", b)
a ist ein <class 'numpy.ndarray'> mit Wert [1. 2. 3.]
b ist ein <class 'numpy.ndarray'> mit Wert [6. 5. 4.]
Die Klasse ndarray repräsentiert ein mehrdimensionales Array. In unserem Fall ein Vektor, also ein Array der Dimension
a.ndim
1
mit Daten vom Typ
a.dtype
dtype('float64')
der Dimension
a.shape
(3,)
2 weitere Funktionen um NumPy-Arrays zu erzeugen sind:
x = np.linspace(0,3,6) # Äquidistantes Punktgitter für [0,3] aus 6 Punkten
x
array([0. , 0.6, 1.2, 1.8, 2.4, 3. ])
x = np.arange(0, 3, 0.5) # Äquidistantes Punktgitter für [0,3) mit Inkrement 0.5
x
array([0. , 0.5, 1. , 1.5, 2. , 2.5])
3.2.2. Elementare Vektoroperationen#
Ein ndarray ist, wie eine Liste oder ein Tupel, ein iterierbares Objekt. Wir können also einfach mit einer for-Schleife über alle Elemente gehen:
val = 0
for e in a:
val += e
print("Die Summe der Einträge von", a, "ist", val)
Die Summe der Einträge von [1. 2. 3.] ist 6.0
Wir können einige elementare Rechenoperationen für mit unseren Vektoren durchführen und es kommt das erwartete Ergebnis raus:
a+b
array([7., 7., 7.])
b-a
array([5., 3., 1.])
a*b
array([ 6., 10., 12.])
a/b
array([0.16666667, 0.4 , 0.75 ])
Wir beobachten, dass die Grundrechenarten einfach elementweise durchgeführt werden. Bei Addition und Subtraktion ist es auch das, was wir aus der Linearen Algebra-Perspektive erwarten würden, aber bei der Multiplikation und Division ist dies nicht klar. Eventuell hätten wir die Berechnung des Skalar- oder Kreuzproduktes erwartet. Dazu später mehr.
Die Rechenoperation +, -, *, / werfen auch einen Fehler, wenn die Größen der beiden Vektoren nicht kompatibel sind:
c = np.array([8,7,6,5])
a+c
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[13], line 2
1 c = np.array([8,7,6,5])
----> 2 a+c
ValueError: operands could not be broadcast together with shapes (3,) (4,)
Untersuchen wir zunächst welche Attribute und Methoden von der Klasse ndarray bereitgestelt werden. Dazu geben wir a.<TAB> bzw. dir(a) ein und bekommen eine lange Liste angezeigt. Testen wir 3 dieser Methoden aus:
print("Der kleinste Eintrag von a ist", a.min(), "bei Index", a.argmin())
print("Das Skalarprodukt <a,b> ist", a.dot(b))
print("Die Summe aller Einträge von a ist", a.sum())
Der kleinste Eintrag von a ist 1.0 bei Index 0
Das Skalarprodukt <a,b> ist 28.0
Die Summe aller Einträge von a ist 6.0
3.2.3. Weitere Rechenoperationen#
Neben diesen elementaren Funktionen, welche die Klasse ndarray direkt bereitstellt, finden wir noch weitere Rechenoperation, welche als freie Funktionen in der Bibliothek numpy implementiert sind. Wir tippen np.<TAB> ein und erhalten eine Liste all dieser Funktionen. Uns fällt auf, dass hier nochmals die Funktionen sqrt, exp, sin, cos definiert sind. Diese sind zwar schon in der math-Bibliothek vorhanden, lassen sich aber nicht auf Objekte vom Typ ndarray anwenden:
import math
math.exp(a)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[15], line 2
1 import math
----> 2 math.exp(a)
TypeError: only length-1 arrays can be converted to Python scalars
Die Exponentialfunktion aus der numpy-Bibliothek hingegen wendet die Exponentialfunktion komponentenweise auf den Vektor an:
np.exp(a)
array([ 2.71828183, 7.3890561 , 20.08553692])
Auch die Funktionen log, sin, cos, tan, abs, …, werden elementweise auf den Vektor angewendet. Genau so werden auch Vergleichsoperationen elementweise angewendet:
x = np.array([0,1,2,3,4])
y = (x <=2)
y
array([ True, True, True, False, False])
Übungsaufgabe
Berechne die Werte von \(\sin(x)\) für \(x\in \{0,\frac\pi6,\frac\pi3,\frac\pi2,\ldots,2\pi\}\).
Neben diesen elementaren Rechenoperationen finden wir im Modul numpy aber auch vektorspezifische Operationen, wie beispielsweise das Skalar- und Kreuzprodukt:
np.dot(a,b)
28.0
np.cross(a,b)
array([-7., 14., -7.])
Bei einer genaueren Betrachtung der Funktionen aus np stellt man fest, dass Funktionen für beispielsweise Vektornormen fehlen. Diese finden wir im Submodul numpy.linalg. Wir tippen np.linalg.<TAB> ein und erhalten wieder eine Liste aller angebotenen Funktionen.
Wir können beispielsweise wie folgt die üblichen Normen berechnen:
print("Euklidische Norm :", np.linalg.norm(a))
print("Maximumnorm :", np.linalg.norm(a, np.Inf))
print("1-Norm :", np.linalg.norm(a, 1))
Euklidische Norm : 3.7416573867739413
Maximumnorm : 3.0
1-Norm : 6.0
Vektoroperationen
Zusammenfassend stellen wir fest, dass die NumPy-Bibliothek sehr viele Rechenoperationen für Vektoren aus der linearen Algebra bereitstellt. Diese sind entweder
Klassenfunktionen von
ndarray(z.B.a.min())Freie Funktionen im Paket
numpy(z.B.np.dot(a,b))Freie Funktionen im Paket
numpy.linalg(z.B.np.linalg.norm(a))
Übungsaufgabe
Schreibe eine Funktion, die den Winkel zweier Vektoren über die Formel
berechnet.
3.2.4. Zugriff auf Vektoreinträge#
Schauen wir uns zuletzt noch an, wie auf einzelne oder mehrere Elemente des Vektors in einem bestimmten Bereich zugegriffen werden kann. Dies geschieht mit dem []-Operator:
a = np.linspace(0,1,11)
a
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
print("Vierter Eintrag:") # Die Zählung des Index' beginnt bei 0
a[3]
Vierter Eintrag:
0.30000000000000004
Wir können auch einen Teil des Arrays extrahieren, indem wir einen Indexbereich i:j angeben:
a[3:6] # Indizes zwischen 3 und 5
array([0.3, 0.4, 0.5])
a[:6] # Indizes bis 5
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5])
a[6:] # Indizes ab 6
array([0.6, 0.7, 0.8, 0.9, 1. ])
Wir können auch Teilvektoren überschreiben, natürlich unter Beachtung der Dimension:
a[1:4] = np.ones((3,)) # Schreibe Einsen in die Einträge 1 bis 3
a[6:8] = np.zeros((2,)) # Schreibe Nullen in die Einträge 6 und 7
a
array([0. , 1. , 1. , 1. , 0.4, 0.5, 0. , 0. , 0.8, 0.9, 1. ])
Nebenbei haben wir hier auch 2 Methoden kennengelernt um Arrays mit Einträgen 0 oder 1 zu initialisieren.
3.3. Arbeiten mit Matrizen#
3.3.1. Erzeugen von Matrizen#
Auch für Matrizen nutzen wir den Datentyp numpy.ndarray. Zur Initialisierung übergeben wir der np.array-Funktion eine Liste bestehend aus Listen bestehend aus Zahlen. NumPy erkennt automatisch, dass wir offensichtlich ein zweidimensionales Array, also eine Matrix, erstellen wollen:
A = np.array([[1.,4.,2.],[2.,0.,3.],[1.,1.,4.]])
A
array([[1., 4., 2.],
[2., 0., 3.],
[1., 1., 4.]])
Um eine Einheitsmatrix zu erzeugen nutzen wir
B = np.eye(3)
B
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
Wir können auch eine leere Matrix erzeugen und die Einträge anschließend direkt eintragen:
C = np.zeros((3,3)) # Das Argument ist ein Tupel und gibt die Größer der Matrix an
C[0,0] = 1
C[1,1] = 2
C[1,2] = 5
C[2,1] = 4
C[2,2] = 1
C
array([[1., 0., 0.],
[0., 2., 5.],
[0., 4., 1.]])
Zu einem gegebenen Vektor kann man auch die zugehörige Diagonalmatrix erzeugen mit
D = np.diag(np.array([1.,2.,3.]))
D
array([[1., 0., 0.],
[0., 2., 0.],
[0., 0., 3.]])
3.3.2. Rechnen mit Matrizen#
Wie schon bei den Vektoren beobachtet werden die Grundrechenarten sowie Funktionen exp, sin, cos, log, etx. aus dem Paket numpy elementweise angewendet:
C = A+B
C
array([[2., 4., 2.],
[2., 1., 3.],
[1., 1., 5.]])
D = A*B
D
array([[1., 0., 0.],
[0., 0., 0.],
[0., 0., 4.]])
E = np.exp(A)
E
array([[ 2.71828183, 54.59815003, 7.3890561 ],
[ 7.3890561 , 1. , 20.08553692],
[ 2.71828183, 2.71828183, 54.59815003]])
Weitere Rechenoperationen werden von der Klasse numpy.ndarray bereitgestellt, wie das Transponieren einer Matrix:
C.transpose()
array([[2., 2., 1.],
[4., 1., 1.],
[2., 3., 5.]])
Die Matrixmultiplikation wird als freie Funktion vom Modul numpy bereitgestellt:
np.matmul(A,B)
array([[1., 4., 2.],
[2., 0., 3.],
[1., 1., 4.]])
oder alternativ mit dem @-Operator
A@B
array([[1., 4., 2.],
[2., 0., 3.],
[1., 1., 4.]])
Auch in der Bibliothek np.linalg finden wir weitere Funktionen aus der linearen Algebra, beispielsweise eine Funktion zur Berechnung der Determinante
np.linalg.det(C)
-22.000000000000004
der Inversen
np.linalg.inv(C)
array([[-0.09090909, 0.81818182, -0.45454545],
[ 0.31818182, -0.36363636, 0.09090909],
[-0.04545455, -0.09090909, 0.27272727]])
und der Eigenwerte- und Vektoren
E,V = np.linalg.eig(C)
print("Eigenwerte :\n", E)
print("Eigenvektoren :\n", V)
Eigenwerte :
[ 6.97413644 -1.33574651 2.36161007]
Eigenvektoren :
[[ 0.63932126 0.77319109 -0.8514755 ]
[ 0.50515119 -0.63379132 -0.29406667]
[ 0.57973321 -0.02200211 0.43418229]]
Die Eigenvektoren stehen spaltenweise in V. Testen wir die Rechnung indem wir \(C\,v - \lambda\,v=0\) überprüfen:
for i in range(3):
res = np.matmul(C, V[:,i]) - E[i]*V[:,i] # C*v-lambda*v
print("Fehler: ", np.linalg.norm(res))
Fehler: 3.66205343881779e-15
Fehler: 1.5852031831917945e-15
Fehler: 2.434909185329438e-15
Auch für die Lösung linearer Gleichungssysteme gibt es eine Funktion aus dem Modul
numpy.linalg:
b = np.array([24.,23.,30.])
x = np.linalg.solve(C, b)
x
array([3., 2., 5.])
Eine einfache Probe ergibt
np.matmul(C,x)-b
array([ 0.00000000e+00, 0.00000000e+00, -3.55271368e-15])
Zugriff auf Matrixeinträge
Auch der Zugriff auf einzelne Einträge bzw. auf Teilmatrizen funktioniert analog zu eindimensionalen NumPy-Arrays:
print("Eintrag C_11 :", C[0,0])
print("Zweite Spalte :", C[:,1])
print("Dritte Zeile :", C[2,:])
print("Letzte (=dritte) Zeile :", C[-1,:])
Eintrag C_11 : 2.0
Zweite Spalte : [4. 1. 1.]
Dritte Zeile : [1. 1. 5.]
Letzte (=dritte) Zeile : [1. 1. 5.]
Übungsaufgabe
Schreibe eine Funktion, welche überprüft, ob eine gegebene Matrix eine Nullzeile besitzt.
Matrizen stapeln
Mit den Befehlen numpy.vstack und numpy.hstack lassen sich Matrizen vertikal bzw. horizontal “zusammenkleben”:
np.vstack([A,B,C]) # Stapelt A, B, C vertikal
array([[1., 4., 2.],
[2., 0., 3.],
[1., 1., 4.],
[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.],
[2., 4., 2.],
[2., 1., 3.],
[1., 1., 5.]])
np.hstack([A,B,C]) # Stapelt A, B, C vertikal
array([[1., 4., 2., 1., 0., 0., 2., 4., 2.],
[2., 0., 3., 0., 1., 0., 2., 1., 3.],
[1., 1., 4., 0., 0., 1., 1., 1., 5.]])
Um einen Vektor an eine Matrix zu kleben müssen wir bei numpy.hstack den Vektor zunächst in eine \(n\times 1\)-Matrix umwandeln wir folgendes Beispiel zeigt:
np.vstack([A,b])
array([[ 1., 4., 2.],
[ 2., 0., 3.],
[ 1., 1., 4.],
[24., 23., 30.]])
np.hstack([A,b.reshape((3,1))])
array([[ 1., 4., 2., 24.],
[ 2., 0., 3., 23.],
[ 1., 1., 4., 30.]])
Mutable oder Immutable?
Zuletzt wollen wir noch überprüfen, wie sich NumPy-Arrays als Funktionsparameter verhalten. Sind diese mutable oder immutable? Ein einfacher Test gibt:
def modify_matrix(A):
A[1,1] = 1.
A = np.zeros((3,3))
modify_matrix(A)
A
array([[0., 0., 0.],
[0., 1., 0.],
[0., 0., 0.]])
Objekte vom Typ numpy.ndarray sind offensichtlich mutable. Man muss an dieser Stelle auch bei einer Zuweisung von Matrizen aufpassen:
B = A
B[1,1] = 2.
A
array([[0., 0., 0.],
[0., 2., 0.],
[0., 0., 0.]])
Wir haben hier den Namen B an das gleiche Objekt wie A gebunden. Nachdem wir B verändert haben, hat sich diese Änderung offensichtlich auch auf die Matrix A ausgewirkt. Es wird also keine Kopie der Matrix angezeigt.
Möchte man tatsächlich eine Kopie einer Matrix erstellen, so nutzt man
B = np.copy(A)
B[1,1] = 3.
A
array([[0., 0., 0.],
[0., 2., 0.],
[0., 0., 0.]])
Wir haben hier B verändert, die Änderung wirkt sich aber nicht auf die Matrix A aus.
Übungsaufgabe
Erstelle die Matrix
welche bei der Finite-Differenzen-Diskretisierung des Randwertproblems \(-y''(t) = f(t)\) für \(t\in(0,1)\) und \(y(0)=y(1)=0\) auftritt.