Lektion 9



In [1]:

    
from sympy import *
init_printing()
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import display

Funktionen als Argumente

"Functions are first class citizens"



In [2]:

    
def komposition(f, g):
    "gibt die Funktion f ∘ g zurück"
    def func(x):
        return f(g(x))
    return func



In [3]:

    
h = komposition(cos, sqrt)
h(pi**2)









    Out[3]:





$$-1$$



In [4]:

    
x = Symbol('x')
h(x)









    Out[4]:





$$\cos{\left (\sqrt{x} \right )}$$

Jordansche Normalform



In [5]:

    
C = Matrix(3, 3, [-4, -2, -3, 5, 3, 3, 5, 2, 4])
C









    Out[5]:





$$\left[\begin{matrix}-4 & -2 & -3\\5 & 3 & 3\\5 & 2 & 4\end{matrix}\right]$$



In [6]:

    
C.eigenvals()









    Out[6]:





$$\left \{ 1 : 3\right \}$$



In [7]:

    
C.eigenvects()









    Out[7]:





$$\left [ \left ( 1, \quad 3, \quad \left [ \left[\begin{matrix}- \frac{2}{5}\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}- \frac{3}{5}\\0\\1\end{matrix}\right]\right ]\right )\right ]$$

Das ist einer zu wenig.



In [8]:

    
T, J = C.diagonalize()  # MatrixError









    



---------------------------------------------------------------------------
MatrixError                               Traceback (most recent call last)
<ipython-input-8-3144553eeb69> in <module>()
----> 1 T, J = C.diagonalize()  # MatrixError

C:\Users\braun\Miniconda3\lib\site-packages\sympy\matrices\matrices.py in diagonalize(self, reals_only, sort, normalize)
   3465         if not self.is_diagonalizable(reals_only, False):
   3466             self._diagonalize_clear_subproducts()
-> 3467             raise MatrixError("Matrix is not diagonalizable")
   3468         else:
   3469             if self._eigenvects is None:

MatrixError: Matrix is not diagonalizable



In [10]:

    
T, J = C.jordan_form()
T, J









    Out[10]:





$$\left ( \left[\begin{matrix}- \frac{38}{3} & \frac{5}{3} & - \frac{1}{7}\\\frac{38}{3} & \frac{2}{3} & - \frac{8}{7}\\\frac{38}{3} & 1 & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 1 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\right )$$



In [11]:

    
T * J * T**(-1) == C









    Out[11]:





True

Rang einer Matrix



In [12]:

    
M = Matrix(3, 3, range(1,10))
M









    Out[12]:





$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$$



In [13]:

    
M.rank()









    Out[13]:





$$2$$



In [14]:

    
x = Symbol('x')
y = Symbol('y')
M = Matrix(3, 2, [2*x+2, 2*y-2, 2*x+2, -2*y+2, y-1, x+1])
M









    Out[14]:





$$\left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & - 2 y + 2\\y - 1 & x + 1\end{matrix}\right]$$



In [15]:

    
M.rank()









    Out[15]:





$$2$$

Glauben wir das?



In [16]:

    
def mult_zeile(A, zeilennr, faktor):
    B = Matrix(A)
    def func(b, j):
        return b*faktor
    B.row_op(zeilennr, func)
    return B

def add_zeile(A, quell_nr, ziel_nr, faktor):
    B = Matrix(A)  
    def func(b, j):
        return b + B[quell_nr,j]*faktor
    B.row_op(ziel_nr, func)
    return B



In [17]:

    
M1 = add_zeile(M, 0, 1, -1)
M1









    Out[17]:





$$\left[\begin{matrix}2 x + 2 & 2 y - 2\\0 & - 4 y + 4\\y - 1 & x + 1\end{matrix}\right]$$



In [18]:

    
M2 = mult_zeile(M1, 0, y-1)   # nur für y != 1
M3 = mult_zeile(M2, 2, x+1)   # nur für x != -1
M3









    Out[18]:





$$\left[\begin{matrix}\left(2 x + 2\right) \left(y - 1\right) & \left(y - 1\right) \left(2 y - 2\right)\\0 & - 4 y + 4\\\left(x + 1\right) \left(y - 1\right) & \left(x + 1\right)^{2}\end{matrix}\right]$$



In [19]:

    
M4 = add_zeile(M3, 0, 2, -Rational(1,2)).expand()
factor(M4)









    Out[19]:





$$\left[\begin{matrix}2 \left(x + 1\right) \left(y - 1\right) & 2 \left(y - 1\right)^{2}\\0 & - 4 \left(y - 1\right)\\0 & \left(x + y\right) \left(x - y + 2\right)\end{matrix}\right]$$

Wenn $x\ne-1$ und $y\ne1$, dann ist der Rang dieser Matrix gleich 2



In [20]:

    
M5 = M.subs(x, -1)
factor(M5)









    Out[20]:





$$\left[\begin{matrix}0 & 2 \left(y - 1\right)\\0 & - 2 \left(y - 1\right)\\y - 1 & 0\end{matrix}\right]$$

Für $y\ne1$ ist der Rang dieser Matrix immer noch 2.



In [21]:

    
M6 = M.subs(y, 1)
factor(M6)









    Out[21]:





$$\left[\begin{matrix}2 \left(x + 1\right) & 0\\2 \left(x + 1\right) & 0\\0 & x + 1\end{matrix}\right]$$

dito für $x\ne-1$



In [22]:

    
M7 = M.subs({x:-1, y:1})
M7









    Out[22]:





$$\left[\begin{matrix}0 & 0\\0 & 0\\0 & 0\end{matrix}\right]$$

Kein Fehler von sympy.

Normen von Vektoren und Matrizen



In [23]:

    
v = Matrix([1,2,3])
v









    Out[23]:





$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$



In [24]:

    
v.norm()









    Out[24]:





$$\sqrt{14}$$



In [27]:

    
sqrt((v.T*v)[0])









    Out[27]:





$$\sqrt{14}$$



In [28]:

    
v.norm(oo)









    Out[28]:





$$3$$



In [30]:

    
v.norm(4)









    Out[30]:





$$\sqrt[4]{2} \sqrt{7}$$



In [32]:

    
C









    Out[32]:





$$\left[\begin{matrix}-4 & -2 & -3\\5 & 3 & 3\\5 & 2 & 4\end{matrix}\right]$$



In [33]:

    
C.norm()









    Out[33]:





$$3 \sqrt{13}$$

Die Frobeniusnorm ist submultiplikativ, also $\left\Vert A B \right\Vert \le \left\Vert A \right\Vert \left\Vert B \right\Vert$. Sie ist aber keine Matrixnorm im Sinne der Analysis II.



In [34]:

    
C.norm(2)









    Out[34]:





$$\sqrt{\sqrt{3363} + 58}$$



In [35]:

    
(C.T*C).eigenvals()









    Out[35]:





$$\left \{ 1 : 1, \quad - \sqrt{3363} + 58 : 1, \quad \sqrt{3363} + 58 : 1\right \}$$



In [36]:

    
C.norm().n(), C.norm(2).n()









    Out[36]:





$$\left ( 10.816653826392, \quad 10.7699293716158\right )$$

Kreuzprodukt



In [37]:

    
w = Matrix([-1,-2,3])
z = v.cross(w)
z









    Out[37]:





$$\left[\begin{matrix}12\\-6\\0\end{matrix}\right]$$



In [38]:

    
v.T*z









    Out[38]:





$$\left[\begin{matrix}0\end{matrix}\right]$$



In [39]:

    
w.T*z









    Out[39]:





$$\left[\begin{matrix}0\end{matrix}\right]$$

Das Kreuzprodukt $v \times w$ steht senkrecht auf $v$ und $w$. Seine euklidische Norm ist gleich dem Flächeninhalt des von $v$ und $w$ aufgespannten Parallelogramms.

Gradient



In [40]:

    
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')
var = [x,y,z]



In [41]:

    
f = exp(x**2 + 2*y**2 + 3*z**2)
f









    Out[41]:





$$e^{x^{2} + 2 y^{2} + 3 z^{2}}$$



In [42]:

    
gr = Matrix([f.diff(t) for t in var])
gr









    Out[42]:





$$\left[\begin{matrix}2 x e^{x^{2} + 2 y^{2} + 3 z^{2}}\\4 y e^{x^{2} + 2 y^{2} + 3 z^{2}}\\6 z e^{x^{2} + 2 y^{2} + 3 z^{2}}\end{matrix}\right]$$

sympy hat eine Funktion gradient, aber nur in 3D. Man beachte aber den nächsten Abschnitt.

Jacobi- und Hessematrix



In [45]:

    
J = Matrix([f]).jacobian(var)
J









    Out[45]:





$$\left[\begin{matrix}2 x e^{x^{2} + 2 y^{2} + 3 z^{2}} & 4 y e^{x^{2} + 2 y^{2} + 3 z^{2}} & 6 z e^{x^{2} + 2 y^{2} + 3 z^{2}}\end{matrix}\right]$$



In [46]:

    
H = factor(gr.jacobian(var))
H









    Out[46]:





$$\left[\begin{matrix}2 \left(2 x^{2} + 1\right) e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 8 x y e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 12 x z e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}}\\8 x y e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 4 \left(4 y^{2} + 1\right) e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 24 y z e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}}\\12 x z e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 24 y z e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}} & 6 \left(6 z^{2} + 1\right) e^{x^{2}} e^{2 y^{2}} e^{3 z^{2}}\end{matrix}\right]$$

Das ist die Hesse Matrix.



In [47]:

    
factor(hessian(f, var)) == H









    Out[47]:





True



In [53]:

    
f = Function('f')
f(*var)









    Out[53]:





$$f{\left (x,y,z \right )}$$



In [54]:

    
hessian(f(x,y,z), var)









    Out[54]:





$$\left[\begin{matrix}\frac{\partial^{2}}{\partial x^{2}}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial x\partial y}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial x\partial z}  f{\left (x,y,z \right )}\\\frac{\partial^{2}}{\partial x\partial y}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial y^{2}}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial y\partial z}  f{\left (x,y,z \right )}\\\frac{\partial^{2}}{\partial x\partial z}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial y\partial z}  f{\left (x,y,z \right )} & \frac{\partial^{2}}{\partial z^{2}}  f{\left (x,y,z \right )}\end{matrix}\right]$$

Definitheit



In [55]:

    
H1 = H.subs({x:1, y:0, z:-1})
H1









    Out[55]:





$$\left[\begin{matrix}6 e^{4} & 0 & - 12 e^{4}\\0 & 4 e^{4} & 0\\- 12 e^{4} & 0 & 42 e^{4}\end{matrix}\right]$$



In [58]:

    
H1.is_positive   # falscher Befehl









    Out[58]:





False



In [59]:

    
H1.det()









    Out[59]:





$$432 e^{12}$$



In [60]:

    
for tmp in H1.eigenvals():
    display((tmp, tmp.n()))









    





$$\left ( - 6 \sqrt{13} e^{4} + 24 e^{4}, \quad 129.217023855597\right )$$






    





$$\left ( 6 \sqrt{13} e^{4} + 24 e^{4}, \quad 2491.49417773533\right )$$






    





$$\left ( 4 e^{4}, \quad 218.392600132577\right )$$



In [61]:

    
ImmutableMatrix(eye(3)).is_positive  # das ist nicht die 
# Antwort auf unsere Frage









    Out[61]:





False

Hurwitz Kriterium



In [62]:

    
for j in range(1,4):
    minor = H1[0:j, 0:j]
    display(minor.det())









    





$$6 e^{4}$$






    





$$24 e^{8}$$






    





$$432 e^{12}$$

Also positiv definit



In [63]:

    
for j in range(1,4):
    minor = -H1[0:j, 0:j]
    display(minor.det())









    





$$- 6 e^{4}$$






    





$$24 e^{8}$$






    





$$- 432 e^{12}$$

Hurwitz-Kriterium

Es sei $M \in \mathbb R^{n\times n}$ eine symmetrische Matrix.

$M$ ist genau dann positiv definit, wenn alle Unterdeterminanten entlang der Hauptdiagonale positiv sind.
$M$ ist genau dann negativ definit, wenn die geraden Unterdeterminanten positiv und die ungeraden negativ sind.
Wenn eine gerade Unterdeterminante negativ ist, dann ist $M$ indefinit.

Extremwerte in mehreren Veränderlichen



In [64]:

    
from mpl_toolkits.mplot3d import Axes3D



In [65]:

    
f = -x**4/2 - x**2*y**2 - y**4/2 + x**3 - 3*x*y**2
fn = lambdify((x,y), f)
f









    Out[65]:





$$- \frac{x^{4}}{2} + x^{3} - x^{2} y^{2} - 3 x y^{2} - \frac{y^{4}}{2}$$



In [66]:

    
xn = np.linspace(-5, 5)
yn = np.linspace(-5, 5)
X, Y = np.meshgrid(xn, yn)
W = fn(X, Y)



In [67]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, W, rstride=1, cstride=1, linewidth=0,
                cmap=plt.cm.coolwarm);



In [68]:

    
glg1 = Eq(f.diff(x), 0)
glg2 = Eq(f.diff(y), 0)
gls = {glg1, glg2}
gls









    Out[68]:





$$\left\{- 2 x^{2} y - 6 x y - 2 y^{3} = 0, - 2 x^{3} + 3 x^{2} - 2 x y^{2} - 3 y^{2} = 0\right\}$$



In [69]:

    
Lsg = solve(gls)
Lsg









    Out[69]:





$$\left [ \left \{ x : - \frac{3}{4}, \quad y : - \frac{3 \sqrt{3}}{4}\right \}, \quad \left \{ x : - \frac{3}{4}, \quad y : \frac{3 \sqrt{3}}{4}\right \}, \quad \left \{ x : 0, \quad y : 0\right \}, \quad \left \{ x : \frac{3}{2}, \quad y : 0\right \}\right ]$$



In [70]:

    
H = hessian(f, [x,y])
H









    Out[70]:





$$\left[\begin{matrix}- 6 x^{2} + 6 x - 2 y^{2} & - 4 x y - 6 y\\- 4 x y - 6 y & - 2 x^{2} - 6 x - 6 y^{2}\end{matrix}\right]$$



In [71]:

    
H1 = H.subs(Lsg[0])
H1









    Out[71]:





$$\left[\begin{matrix}- \frac{45}{4} & \frac{9 \sqrt{3}}{4}\\\frac{9 \sqrt{3}}{4} & - \frac{27}{4}\end{matrix}\right]$$



In [72]:

    
H1[0,0]









    Out[72]:





$$- \frac{45}{4}$$



In [73]:

    
H1.det()









    Out[73]:





$$\frac{243}{4}$$

negativ definit



In [74]:

    
H2 = H.subs(Lsg[1])
H2









    Out[74]:





$$\left[\begin{matrix}- \frac{45}{4} & - \frac{9 \sqrt{3}}{4}\\- \frac{9 \sqrt{3}}{4} & - \frac{27}{4}\end{matrix}\right]$$



In [75]:

    
H2.det()









    Out[75]:





$$\frac{243}{4}$$



In [76]:

    
H3 = H.subs(Lsg[2])
H3









    Out[76]:





$$\left[\begin{matrix}0 & 0\\0 & 0\end{matrix}\right]$$



In [77]:

    
H4 = H.subs(Lsg[3])
H4









    Out[77]:





$$\left[\begin{matrix}- \frac{9}{2} & 0\\0 & - \frac{27}{2}\end{matrix}\right]$$



In [78]:

    
for l in Lsg:
    y0 = f.subs(l)
    display((y0, y0.n()))









    





$$\left ( \frac{27}{32}, \quad 0.84375\right )$$






    





$$\left ( \frac{27}{32}, \quad 0.84375\right )$$






    





$$\left ( 0, \quad 0\right )$$






    





$$\left ( \frac{27}{32}, \quad 0.84375\right )$$



In [79]:

    
from matplotlib.colors import Normalize
norm = Normalize(-.5, 1)



In [80]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, W, rstride=1, cstride=1, linewidth=0,
                cmap=plt.cm.coolwarm, norm=norm);



In [81]:

    
for l in Lsg:
    display((x.subs(l).n(), y.subs(l).n()))









    





$$\left ( -0.75, \quad -1.29903810567666\right )$$






    





$$\left ( -0.75, \quad 1.29903810567666\right )$$






    





$$\left ( 0, \quad 0\right )$$






    





$$\left ( 1.5, \quad 0\right )$$



In [82]:

    
xn = np.linspace(-1.5, 2, 100)
yn = np.linspace(-1.8, 1.8, 100)
X, Y = np.meshgrid(xn, yn)
W = np.maximum(fn(X, Y), -.1)



In [83]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, W, rstride=1, cstride=1, linewidth=0,
                cmap=plt.cm.coolwarm);



In [84]:

    
f.subs(y, 0)









    Out[84]:





$$- \frac{x^{4}}{2} + x^{3}$$

Also Sattel in $(0,0)$.



In [ ]: