Computergestützte Mathematik zur Analysis¶

Vorlesung vom 15.12.2022

© 2022 Prof. Dr. Rüdiger W. Braun

In [1]:
from sympy import *
init_printing()

Spezielle Matrizen¶

In [2]:
M = eye(5)
M
Out[2]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 1\end{matrix}\right]$
In [3]:
ones(3,4)
Out[3]:
$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1\\1 & 1 & 1 & 1\\1 & 1 & 1 & 1\end{matrix}\right]$
In [4]:
zeros(2,1)
Out[4]:
$\displaystyle \left[\begin{matrix}0\\0\end{matrix}\right]$
In [5]:
diag(1,2,3)
Out[5]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$
In [6]:
l = [1,2,3]
In [7]:
diag(l)
Out[7]:
$\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]$
In [8]:
diag(*l)
Out[8]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$
In [9]:
print(l)
print(*l)   # *l packt die ELemente von l aus

[1, 2, 3]
1 2 3

Operationen¶

In [11]:
A = Matrix(3, 3, list(range(1,10)))
A
Out[11]:
$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$
In [12]:
A.det()
Out[12]:
$\displaystyle 0$
In [13]:
C = A + eye(3)
C
Out[13]:
$\displaystyle \left[\begin{matrix}2 & 2 & 3\\4 & 6 & 6\\7 & 8 & 10\end{matrix}\right]$
In [14]:
C.det()
Out[14]:
$\displaystyle -2$
In [15]:
C1 = C**(-1)
C1
Out[15]:
$\displaystyle \left[\begin{matrix}-6 & -2 & 3\\-1 & \frac{1}{2} & 0\\5 & 1 & -2\end{matrix}\right]$
In [16]:
C1 == C.inv()
Out[16]:
True
In [17]:
C * C1
Out[17]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

Inspektion¶

In [18]:
A.shape
Out[18]:
$\displaystyle \left( 3, \ 3\right)$
In [19]:
v = Matrix([1,2,3])
v
Out[19]:
$\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]$
In [20]:
v.shape
Out[20]:
$\displaystyle \left( 3, \ 1\right)$

Es gibt keine Vektoren, nur $n \times 1$-Matrizen

In [21]:
type(v)
Out[21]:
sympy.matrices.dense.MutableDenseMatrix

Slicing¶

In [22]:
B = Matrix(5,6, range(30))
B
Out[22]:
$\displaystyle \left[\begin{matrix}0 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right]$
In [23]:
B[2:4, 3:4]
Out[23]:
$\displaystyle \left[\begin{matrix}15\\21\end{matrix}\right]$
In [25]:
B[:4, 3:]
Out[25]:
$\displaystyle \left[\begin{matrix}3 & 4 & 5\\9 & 10 & 11\\15 & 16 & 17\\21 & 22 & 23\end{matrix}\right]$
In [26]:
B[:, :2]
Out[26]:
$\displaystyle \left[\begin{matrix}0 & 1\\6 & 7\\12 & 13\\18 & 19\\24 & 25\end{matrix}\right]$

Kopien¶

In [27]:
C = B
In [28]:
C[0,0] = 121
B, C
Out[28]:
$\displaystyle \left( \left[\begin{matrix}121 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right], \ \left[\begin{matrix}121 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right]\right)$
In [29]:
C = B.copy()
C[0,1] = -1000
B, C
Out[29]:
$\displaystyle \left( \left[\begin{matrix}121 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right], \ \left[\begin{matrix}121 & -1000 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right]\right)$
In [30]:
d1 = {'b':B}  # dictionary
d2 = d1.copy()
In [31]:
d1['b'][0,0] = 42
In [32]:
d1['b'], d2['b']
Out[32]:
$\displaystyle \left( \left[\begin{matrix}42 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right], \ \left[\begin{matrix}42 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right]\right)$
In [33]:
import copy
d3= copy.deepcopy(d1)
d1['b'][0,0] = 12
d1['b'], d3['b']
Out[33]:
$\displaystyle \left( \left[\begin{matrix}12 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right], \ \left[\begin{matrix}42 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\\24 & 25 & 26 & 27 & 28 & 29\end{matrix}\right]\right)$

Manipulation von Matrizen¶

In [34]:
A
Out[34]:
$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$
In [35]:
A.T
Out[35]:
$\displaystyle \left[\begin{matrix}1 & 4 & 7\\2 & 5 & 8\\3 & 6 & 9\end{matrix}\right]$
In [36]:
Matrix.hstack(eye(3), eye(3))   
#  eine  class method
Out[36]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 1 & 0 & 0\\0 & 1 & 0 & 0 & 1 & 0\\0 & 0 & 1 & 0 & 0 & 1\end{matrix}\right]$
In [37]:
Matrix.vstack(eye(3), eye(3))
Out[37]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\\1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$
In [40]:
B = Matrix(4, 6, list(range(24)))
B
Out[40]:
$\displaystyle \left[\begin{matrix}0 & 1 & 2 & 3 & 4 & 5\\6 & 7 & 8 & 9 & 10 & 11\\12 & 13 & 14 & 15 & 16 & 17\\18 & 19 & 20 & 21 & 22 & 23\end{matrix}\right]$
In [42]:
B.reshape(1, 24)
Out[42]:
$\displaystyle \left[\begin{array}{cccccccccccccccccccccccc}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23\end{array}\right]$
In [43]:
flatten(B)
Out[43]:
$\displaystyle \left[ 0, \ 1, \ 2, \ 3, \ 4, \ 5, \ 6, \ 7, \ 8, \ 9, \ 10, \ 11, \ 12, \ 13, \ 14, \ 15, \ 16, \ 17, \ 18, \ 19, \ 20, \ 21, \ 22, \ 23\right]$
In [45]:
type(B.reshape(1,24))
Out[45]:
sympy.matrices.dense.MutableDenseMatrix
In [46]:
type(flatten(B))
Out[46]:
list

Hilbert-Matrizen¶

als Beispiel für schlecht konditionierte Matrizen

In [47]:
def hilbert(i,j):
    return 1/(1+i+j)
In [48]:
H = Matrix(5, 5, hilbert)
H
Out[48]:
$\displaystyle \left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\\frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8}\\\frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9}\end{matrix}\right]$
In [49]:
H.det()
Out[49]:
$\displaystyle \frac{1}{266716800000}$
In [50]:
H.inv()
Out[50]:
$\displaystyle \left[\begin{matrix}25 & -300 & 1050 & -1400 & 630\\-300 & 4800 & -18900 & 26880 & -12600\\1050 & -18900 & 79380 & -117600 & 56700\\-1400 & 26880 & -117600 & 179200 & -88200\\630 & -12600 & 56700 & -88200 & 44100\end{matrix}\right]$
In [51]:
N = 12
In [52]:
H = Matrix(N, N, hilbert)
H.det()
Out[52]:
$\displaystyle \frac{1}{379106579436304517151885479034796391880188687864118464104324304732160000000000}$
In [53]:
H1 = H**(-1)
H1[N-5:, N-5:]
Out[53]:
$\displaystyle \left[\begin{matrix}2654320456719360 & -3110531785218000 & 2276990587872000 & -946216088737920 & 170392979877120\\-3110531785218000 & 3659449159080000 & -2688113888460000 & 1120519052452800 & -202341663604080\\2276990587872000 & -2688113888460000 & 1980715496760000 & -827939077645680 & 149882713780800\\-946216088737920 & 1120519052452800 & -827939077645680 & 346945899203904 & -62950739787936\\170392979877120 & -202341663604080 & 149882713780800 & -62950739787936 & 11445589052352\end{matrix}\right]$

Vergleich mit Numerik¶

In [56]:
import numpy as np
In [57]:
Hn = np.empty((N,N))
for i in range(N):
    for j in range(N):
        Hn[i,j] = 1/(1+i+j)
np.linalg.det(Hn)
Out[57]:
$\displaystyle 2.76978138495454 \cdot 10^{-78}$
In [58]:
H.det().n()
Out[58]:
$\displaystyle 2.63778065125355 \cdot 10^{-78}$
In [60]:
Hn1 = Hn**(-1)
Hn1[N-5:, N-5:]
Out[60]:
array([[15., 16., 17., 18., 19.],
       [16., 17., 18., 19., 20.],
       [17., 18., 19., 20., 21.],
       [18., 19., 20., 21., 22.],
       [19., 20., 21., 22., 23.]])

Hn**(-1) hat die Kehrwerte aller Matrixeinträge bestimmt

In [61]:
Hn1 = np.linalg.inv(Hn)
Hn1[N-3:, N-3:]
Out[61]:
array([[ 1.88806550e+15, -7.89466867e+14,  1.42959120e+14],
       [-7.89478678e+14,  3.30975528e+14, -6.00766591e+13],
       [ 1.42963152e+14, -6.00774507e+13,  1.09285030e+13]])

zum Vergleich

In [62]:
H1[N-3:, N-3:]
Out[62]:
$\displaystyle \left[\begin{matrix}1980715496760000 & -827939077645680 & 149882713780800\\-827939077645680 & 346945899203904 & -62950739787936\\149882713780800 & -62950739787936 & 11445589052352\end{matrix}\right]$
In [63]:
%%timeit
H**(-1)

152 ms ± 946 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [64]:
%%timeit
Hn1 = np.linalg.inv(Hn)

6.33 µs ± 50.3 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Eine Python-Bibliothek, die numerische Operationen in beliebige wählbarer Genauigkeit durchführen kann, ist mpmath

Rang einer Matrix¶

In [65]:
A
Out[65]:
$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$
In [66]:
A.rank()
Out[66]:
$\displaystyle 2$
In [67]:
x = S('x')
y = S('y')
In [68]:
M = Matrix(3, 2, [2*x+2, 2*y-2, 2*x+2, -2*y+2, y-1, x+1])
M
Out[68]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & 2 - 2 y\\y - 1 & x + 1\end{matrix}\right]$
In [69]:
M.rank()
Out[69]:
$\displaystyle 2$

Glauben wir das für alle Wahlen von $x$ und $y$?

In [70]:
M.rref(pivots=False)  # Zeilenstufenform (engl. reduced row echelon form)
Out[70]:
$\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\\0 & 0\end{matrix}\right]$
In [71]:
M
Out[71]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & 2 - 2 y\\y - 1 & x + 1\end{matrix}\right]$
In [72]:
M1 = M.elementary_row_op('n->kn', row=2, k=2*x+2)
M1
Out[72]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & 2 - 2 y\\\left(2 x + 2\right) \left(y - 1\right) & \left(x + 1\right) \left(2 x + 2\right)\end{matrix}\right]$

Das darf ich aber nur, wenn $2x+2\ne0$.

In [73]:
M2 = M1.elementary_row_op('n->n+km', row1=2, row2=0, k=1-y).expand()
M2
Out[73]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & 2 - 2 y\\0 & 2 x^{2} + 4 x - 2 y^{2} + 4 y\end{matrix}\right]$

Das darf ich immer.

In [74]:
M2
Out[74]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\2 x + 2 & 2 - 2 y\\0 & 2 x^{2} + 4 x - 2 y^{2} + 4 y\end{matrix}\right]$
In [75]:
M3 = M2.elementary_row_op('n->n+km', row1=1, row2=0, k=-1)
M3
Out[75]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\0 & 4 - 4 y\\0 & 2 x^{2} + 4 x - 2 y^{2} + 4 y\end{matrix}\right]$
In [76]:
M4 = M3.elementary_row_op('n->n+km', row1=2, row2=1, k=-M3[2,1]/M3[1,1])
M4
Out[76]:
$\displaystyle \left[\begin{matrix}2 x + 2 & 2 y - 2\\0 & 4 - 4 y\\0 & 0\end{matrix}\right]$

Das darf ich nur für $4-4y\ne0$, weil ich durch diesen Wert geteilt habe

Bis jetzt gesehen:

Für $x\ne-1$ und $y\ne1$ ist der Rang gleich $2$.

In [77]:
M.subs(x, -1)
Out[77]:
$\displaystyle \left[\begin{matrix}0 & 2 y - 2\\0 & 2 - 2 y\\y - 1 & 0\end{matrix}\right]$
In [78]:
M.subs({x:-1, y:1})
Out[78]:
$\displaystyle \left[\begin{matrix}0 & 0\\0 & 0\\0 & 0\end{matrix}\right]$
$$ \text{Rang}(M) = \begin{cases} 0, & x = -1 \wedge y = 1, \\ 2, & \text{in allen anderen Fällen.} \end{cases} $$

sympy rechnet im Körper $\mathbb R(x,y)$ der rationalen Polynome