diff --git a/lu.py b/lu.py index a2dd1821147e2e2cfab1e9c4e620c7aa5f17a98d..7ceea1fa2046f077e6888d48bc3789031c3c5a8b 100644 --- a/lu.py +++ b/lu.py @@ -82,7 +82,21 @@ def lu_decompose(A): P = np.eye(n) L = np.zeros((n, n)) U = A.copy() - pass # TODO + # TODO + L = np.eye(n) + for i in range(n): + pivot_id = pivot_index(U, i) + if i != pivot_id: + U = swap_rows(U, i, pivot_id) + # swap columns of P + P[:, [i, pivot_id]] = P[:, [pivot_id, i]] + # swap elts in indexs of L + L[[i, pivot_id], 0:i] = L[[pivot_id, i], 0:i] + for j in range(i + 1, n): + # set scalar in L + L[j][i] = U[j][i] / U[i][i] + # sub Zeilen + subtract_scaled(U, j, i, L[j][i]) return (P, L, U) @@ -96,7 +110,15 @@ def forward_substitute(L, b): Returns: np.ndarray: Vector ``y``, such that (*) holds. """ - pass # TODO + # TODO + n, m = np.shape(L) + y = np.zeros((n, 1)) + y[0] = b[0] / L[0][0] + for i in range(1, n): + for j in range(0, i): + y[i] += L[i][j] * y[j] + y[i] = b[i] - y[i] + return y def backward_substitute(R, y): @@ -109,7 +131,16 @@ def backward_substitute(R, y): Returns: np.ndarray: Vector ``x``, such that (**) holds. """ - pass # TODO + # TODO + n, m = np.shape(R) + x = np.zeros(np.shape(y)) + x[n - 1] = y[n - 1] / R[n - 1][n - 1] + for i in range(n - 2, -1, -1): + for j in range(n - 1, i, -1): + x[i] += R[i][j] * x[j] + x[i] = (y[i] - x[i]) / R[i][i] + return x + def linsolve(A, *bs): @@ -123,7 +154,18 @@ def linsolve(A, *bs): Returns: Tuple[np.ndarray, ...]: ``(x_{0}, ..., x_{n - 1})`` as mentioned above. """ - pass # TODO + # TODO + result = [] + n, m = np.shape(A) + for b in bs: + P, L, U = lu_decompose(A) + if (not np.array_equal(P, np.eye(n))): + transpose = np.transpose(P) + b = np.dot(transpose, b) + y = forward_substitute(L, b) + x = backward_substitute(U, y) + result.append(x) + return tuple(result) def main(): diff --git a/qr.py b/qr.py index 64fcd9c988f2a6cc8e1a3fd434f3e853da8bd91a..ebe08b63c8c16c26dcc663c17c304953661ad963 100644 --- a/qr.py +++ b/qr.py @@ -15,7 +15,16 @@ def givens_rotation(A, i, j): Returns: np.ndarray: Rotation Matrix ``J`` as above. """ - pass # TODO + # TODO + c = abs(A[j][j]) / ((A[i][j] ** 2 + A[j][j] ** 2) ** 0.5) + s = abs(A[i][j]) / ((A[i][j] ** 2 + A[j][j] ** 2) ** 0.5) + n, m = np.shape(A) + J = np.eye(n) + J[j][j] = c + J[j][i] = s + J[i][j] = -s + J[i][i] = c + return J def qr_decompose(A): @@ -27,7 +36,16 @@ def qr_decompose(A): Returns: Tuple[np.ndarray,np.ndarray]: QR-decomposition in order ``(Q, R)``. """ - pass # TODO + # TODO + n, m = np.shape(A) + R = A.copy() + J = np.eye(n) + for j in range(0, n - 1): + for i in range(j + 1, n): + J_new = givens_rotation(R, i, j) + J = np.dot(J_new, J) + R = np.dot(J_new, R) + return (np.transpose(J), R) def backward_substitute(R, y): @@ -40,7 +58,15 @@ def backward_substitute(R, y): Returns: np.ndarray: Vector ``x``, such that (**) holds. """ - pass # TODO + # TODO + n, m = np.shape(R) + x = np.zeros(np.shape(y)) + x[n - 1] = y[n - 1] / R[n - 1][n - 1] + for i in range(n - 2, -1, -1): + for j in range(n - 1, i, -1): + x[i] += R[i][j] * x[j] + x[i] = (y[i] - x[i]) / R[i][i] + return x def linsolve(A, *bs): @@ -54,7 +80,17 @@ def linsolve(A, *bs): Returns: Tuple[np.ndarray, ...]: ``(x_{0}, ..., x_{n - 1})`` as mentioned above. """ - pass # TODO + # TODO + result = [] + n, m = np.shape(A) + for b in bs: + Q, R = qr_decompose(A) + if (not np.array_equal(Q, np.eye(n))): + transpose = np.transpose(Q) + b = np.dot(transpose, b) + x = backward_substitute(R, b) + result.append(x) + return tuple(result) def main():