Gauss-Elimination¶
Berechne das lineare Gleichungssystem
\[\begin{split}\begin{pmatrix}
2 & -1 & -3 & 3\\
4 & 0 & -3 & 1\\
6 & 1 & -1 & 6\\
-2& -5 & 4 & 1\end{pmatrix}\cdot x = \begin{pmatrix}1\\ -8\\ -16\\ -12\end{pmatrix}\end{split}\]
[1]:
import numpy as np
[2]:
A = np.array([[2, -1, -3, 3],
[4, 0, -3, 1],
[6, 1, -1, 6],
[-2, -5, 4, 1]])
Aorig = np.array(A)
A
[2]:
array([[ 2, -1, -3, 3],
[ 4, 0, -3, 1],
[ 6, 1, -1, 6],
[-2, -5, 4, 1]])
[3]:
b = np.array([1, -8, -16, -12])
b
[3]:
array([ 1, -8, -16, -12])
In der Matrix \(L\) werden die Faktoren für die Elimination pro Zeile bzw. Spalte an den analogen Positionen gespeichert, daher unterhalb der Diagonalen.
[4]:
L = np.eye(4)
L
[4]:
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]])
1. Spalte¶
2. Zeile¶
[5]:
spalte = 0
[6]:
zeile = 1
[7]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[7]:
array([[1., 0., 0., 0.],
[2., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]])
Elimination in der Zeile:
[8]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[9]:
A
[9]:
array([[ 2, -1, -3, 3],
[ 0, 2, 3, -5],
[ 6, 1, -1, 6],
[-2, -5, 4, 1]])
3. Zeile¶
[10]:
zeile = 2
[11]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[11]:
array([[1., 0., 0., 0.],
[2., 1., 0., 0.],
[3., 0., 1., 0.],
[0., 0., 0., 1.]])
Elimination in der Zeile:
[12]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[13]:
A
[13]:
array([[ 2, -1, -3, 3],
[ 0, 2, 3, -5],
[ 0, 4, 8, -3],
[-2, -5, 4, 1]])
4. Zeile¶
[14]:
zeile = 3
[15]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[15]:
array([[ 1., 0., 0., 0.],
[ 2., 1., 0., 0.],
[ 3., 0., 1., 0.],
[-1., 0., 0., 1.]])
Elimination in der Zeile:
[16]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[17]:
A
[17]:
array([[ 2, -1, -3, 3],
[ 0, 2, 3, -5],
[ 0, 4, 8, -3],
[ 0, -6, 1, 4]])
2. Spalte¶
3. Zeile¶
[18]:
spalte = 1
[19]:
zeile = 2
[20]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[20]:
array([[ 1., 0., 0., 0.],
[ 2., 1., 0., 0.],
[ 3., 2., 1., 0.],
[-1., 0., 0., 1.]])
Elimination in der Zeile:
[21]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[22]:
A
[22]:
array([[ 2, -1, -3, 3],
[ 0, 2, 3, -5],
[ 0, 0, 2, 7],
[ 0, -6, 1, 4]])
4. Zeile¶
[23]:
zeile = 3
[24]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[24]:
array([[ 1., 0., 0., 0.],
[ 2., 1., 0., 0.],
[ 3., 2., 1., 0.],
[-1., -3., 0., 1.]])
Elimination in der Zeile:
[25]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[26]:
A
[26]:
array([[ 2, -1, -3, 3],
[ 0, 2, 3, -5],
[ 0, 0, 2, 7],
[ 0, 0, 10, -11]])
3. Spalte¶
4. Zeile¶
[27]:
spalte = 2;
[28]:
zeile = 3;
[29]:
L[zeile, spalte] = A[zeile, spalte]/A[spalte, spalte]
L
[29]:
array([[ 1., 0., 0., 0.],
[ 2., 1., 0., 0.],
[ 3., 2., 1., 0.],
[-1., -3., 5., 1.]])
Elimination in der Zeile:
[30]:
for j in range(spalte,4):
A[zeile, j] -= L[zeile, spalte] * A[spalte, j]
[31]:
R = A
[32]:
L
[32]:
array([[ 1., 0., 0., 0.],
[ 2., 1., 0., 0.],
[ 3., 2., 1., 0.],
[-1., -3., 5., 1.]])
Faktorisierung der Matrix¶
\[A = L\cdot R\]
[33]:
L@R
[33]:
array([[ 2., -1., -3., 3.],
[ 4., 0., -3., 1.],
[ 6., 1., -1., 6.],
[-2., -5., 4., 1.]])
[34]:
Aorig
[34]:
array([[ 2, -1, -3, 3],
[ 4, 0, -3, 1],
[ 6, 1, -1, 6],
[-2, -5, 4, 1]])
[35]:
L@R-Aorig
[35]:
array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])
Lösen des Gleichungssystems¶
\[A \cdot x = L\cdot \underbrace{R \cdot x}_{= z} = b\]
Im ersten Schritt wird das System
\[L\cdot z = b\]durch Vorwärtseinsetzen gelöst.
Im zweiten Schritt wird das System
\[R\cdot x = z\]durch Rückwärtseinsetzen gelöst.
Vorwärtseinsetzen¶
\[z_j = b_j - \sum_{i=1}^{j-1} L_{j,i}\cdot z_i,\quad j = 1, \ldots, N\]
[36]:
z = np.zeros(4)
Da \(L\) normiert ist, folgt für die erste Zeile
[37]:
z[0] = b[0]
Für den Resten gilt
[38]:
for zeile in range(1,4):
z[zeile] = b[zeile]-L[zeile,:zeile]@z[:zeile]
[39]:
z
[39]:
array([ 1., -10., 1., -46.])
Rückwärtseinsetzen¶
\[x_j = \frac{z_j - \sum_{i=j+1}^{N} R_{j,i}\cdot x_i}{R_{j,j}},\quad j = N, \ldots, 1\]
[40]:
x = np.zeros(4)
[41]:
x[3] = z[3]/R[3,3]
[42]:
for zeile in range(2,-1,-1):
x[zeile] = (z[zeile]-R[zeile,zeile+1:]@x[zeile+1:])/R[zeile,zeile]
[43]:
x
[43]:
array([-4.5, 2. , -3. , 1. ])
Residuum der Lösung:
[44]:
Aorig@x-b
[44]:
array([0., 0., 0., 0.])