7.2.1. Beispiele zu quadratischen Formen#
import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import solve_triangular, eig, eigvals, norm, qr, svd
from sympy import symbols, lambdify
Mesh Grid für die Visualisierung
xp = np.linspace(-2,2,30)
Xp,Yp = np.meshgrid(xp,xp)
Allgemeine quadratische Form
wobei \(a\in\mathbb{R}^{n\times n}\), \(b\in\mathbb{R}^n\) und \(c\in\mathbb{R}\) sei.
def getQuadForm(a,b,c, lam = True):
x,y = symbols('x,y')
u = np.array([x,y])
expr=u@a@u+b@u+c
if lam:
return lambdify((x,y),expr)
else:
return expr
Elliptisch#
Die Matrix \(a\) muss positiv definit sein:
a = np.array([[2,1],[0,1]])
b = np.array([0,0])
c = 0
getQuadForm(a,b,c, lam = False)
plt.contour(Xp,Yp,getQuadForm(a,b,c)(Xp,Yp))
plt.gca().set_aspect(1)
plt.grid()
plt.show()
a = np.array([[-2,1],[0,-1]])
b = np.array([0,0])
c = 0
getQuadForm(a,b,c, lam = False)
eigvals(a)
array([-2.+0.j, -1.+0.j])
plt.contour(Xp,Yp,getQuadForm(a,b,c)(Xp,Yp))
plt.gca().set_aspect(1)
plt.grid()
plt.show()
Parabolisch#
Ein kleiner Ausrutscher in die lineare Algebra: Berechnung des Kern einer Matrix#
Der Kern einer linearen Abbildung \(A: \mathbb{R}^n \to \mathbb{R}^m\) ist definiert durch die Menge der Vektoren, welche die lineare Abbildung auf den Nullvektor abbildet. Daher
Mit Hilfe der Eigenwertzerlegung#
Eine numerische Approximation des Kern einer (quadratischen) Matrix können wir mit Hilfe der zum Eigenwert \(\lambda = 0\) zugehörigen Eigenvektoren bestimmen. Das funktioniert solange die Vielfachheit eins ist. Im Fall einer Vielfachheit, ist der Eigenvektor komplexkonjugiert, womit wir mit dem Ansatz nicht alle Basisvektoren für den Nullraum erhalten. In dem Fall müsste man die Hauptvektoren bestimmen.
Als erstes erstellen wir eine 5x5 Matrix mit Rang 3:
a = np.random.randint(-9,9,size=(5,3))
# build quadratic matrix with lower rank
a = np.vstack((a.T,a[:,0]-a[:,1],a[:,2]-a[:,1])).T
a
array([[ 1, -1, 0, 2, 1],
[-1, 1, 7, -2, 6],
[ 3, -2, -5, 5, -3],
[-9, -4, -9, -5, -5],
[ 5, -9, 2, 14, 11]])
# Beispiel einer Matrix mit komplex konjugierten Eigenvektoren für den Eigenwert 0.
#a = np.array([[ 3, -7, 6, 10, 13],
# [ 8, 1, 8, 7, 7],
# [ 7, -8, -5, 15, 3],
# [ 2, -7, -8, 9, -1],
# [-5, -6, -8, 1, -2]])
Für die Eigenwerte bzw. Eigenvektoren zum Eigenwert 0 folgt
ew,ev = eig(a)
ns1=ev[:,(np.abs(ew)<1e-13)] # wir sind für das weitere nur am Realanteil interessiert
ns1
array([[ 0.52459553+0.j, -0.00069488+0.j],
[-0.64757376+0.j, -0.57688664+0.j],
[ 0.12297822+0.j, 0.57758152+0.j],
[-0.52459553+0.j, 0.00069488+0.j],
[-0.12297822+0.j, -0.57758152+0.j]])
a@ns1
array([[-6.93889390e-17+0.j, -6.66133815e-16+0.j],
[-2.52575738e-15+0.j, 6.88338275e-15+0.j],
[ 1.31838984e-15+0.j, -6.55031585e-15+0.j],
[ 1.45716772e-15+0.j, 1.66533454e-15+0.j],
[-3.10862447e-15+0.j, -4.44089210e-16+0.j]])
[norm(a@ni,np.inf) for ni in ns1.T]
[np.float64(2.55351295663786e-15), np.float64(7.105427357601002e-15)]
Mit Hilfe der Singulärwertzerlegung#
Ein analoger Zugang, jedoch mit der Singulärwertzerlegung liefert immer eine reellwertige Lösung und den gesammten Nullraum:
atol=1e-13 # absolute Toleranz
rtol=0 # relative Toleranz
u, s, vh = svd(a) # Singulärwerte sind
tol = max(atol, rtol * s[0])
nnz = (s >= tol).sum()
ns2 = vh[nnz:].conj().T
ns2
array([[ 0.61227801, 0.01075368],
[-0.39804667, -0.58443036],
[-0.21423133, 0.57367668],
[-0.61227801, -0.01075368],
[ 0.21423133, -0.57367668]])
a@ns2
array([[2.49800181e-16, 0.00000000e+00],
[6.10622664e-16, 0.00000000e+00],
[1.38777878e-16, 0.00000000e+00],
[1.41553436e-15, 4.44089210e-16],
[2.66453526e-15, 0.00000000e+00]])
[norm(a@ni,np.inf) for ni in ns2.T]
[np.float64(2.6645352591003757e-15), np.float64(2.220446049250313e-16)]
Wir zeigen, dass die Vektoren jeweilen durch Linearkombinationen der anderen darstellbar sind. Womit gezeigt ist, dass die Beschreibungen des Nullraums identisch sind.
q,r = qr(ns1,mode='economic')
r
array([[-1. +0.j, -0.51490749+0.j],
[ 0. +0.j, 0.85724575+0.j]])
sol=solve_triangular(r,q.T@ns2)
sol
array([[ 1.16632276-0.j, 0.02180849-0.j],
[-0.61924355+0.j, 0.98859588+0.j]])
ns1@sol-ns2
array([[ 6.66133815e-16+0.j, -5.79397641e-16+0.j],
[ 4.99600361e-16+0.j, -2.22044605e-16+0.j],
[-4.16333634e-16+0.j, 4.44089210e-16+0.j],
[ 2.22044605e-16+0.j, -3.22658567e-16+0.j],
[-6.93889390e-16+0.j, 4.44089210e-16+0.j]])
Der Vorteil der Singulärwertzerlegung sieht man im Resultat der QR-Zerlegung: die beiden Vektoren sind orthogonal. Die \(R\) Matrix ist diagonal.
q,r = qr(ns2,mode='economic')
r
array([[-1.00000000e+00, 3.51090887e-17],
[ 0.00000000e+00, 1.00000000e+00]])
sol=solve_triangular(r,q.T@ns1)
sol
array([[ 0.84746962-0.j, -0.01869523-0.j],
[ 0.5308439 +0.j, 0.99982523+0.j]])
ns2@sol-ns1
array([[-3.33066907e-16+0.j, 5.19549681e-16+0.j],
[ 0.00000000e+00+0.j, 3.33066907e-16+0.j],
[ 1.38777878e-17+0.j, -5.55111512e-16+0.j],
[ 1.11022302e-16+0.j, 3.96817995e-16+0.j],
[ 4.57966998e-16+0.j, -3.33066907e-16+0.j]])
Die beiden Beschreibungen sind bis auf numerische Rundung identisch.
Zurück zur quadratischen Form#
Betrachten wir das erwähnte Beispiel: Der Nullraum ist in dem Fall offensichtlich gegeben durch \((0,1)^T\) und entsprechend kompatibel mit dem \(b\) Vektor.
a = np.array([[1,0],[0,0]])
b = np.array([0,1])
c = 0
getQuadForm(a,b,c, lam = False)
plt.contour(Xp,Yp,getQuadForm(a,b,c)(Xp,Yp))
plt.gca().set_aspect(1)
plt.grid()
plt.show()
Hyperbolisch#
Wir betrachten das Beispiel:
a = np.array([[1,0],[0,-1]])
b = np.array([0,0])
c = 0
getQuadForm(a,b,c, lam = False)
plt.contour(Xp,Yp,getQuadForm(a,b,c)(Xp,Yp))
plt.gca().set_aspect(1)
plt.grid()
plt.show()
