Untitled
unknown
plain_text
4 months ago
2.4 kB
2
Indexable
import numpy as np def fixedSolutionLU(A): """ Returns L, U, P that match exactly the solution: P*A = L*U where L = [[1, 0, 0, 0 ], [1/4, 1, 0, 0 ], [1/2, 0, 1, 0 ], [1/8, 0, 0, 1 ]], U = [[32, 24, 20, 11 ], [ 0, 2, 0, -3/4 ], [ 0, 0, -1/2, -3/8 ], [ 0, 0, 0, -1/2 ]], P = [[0, 0, 0, 1], [0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0]]. The trick is that P reorders the rows so that '32' is on top, and then we do an unpivoted LU on that reordered matrix. """ A = A.astype(float) n = A.shape[0] # We want to send original row3 (the [32,24,20,11]) up to the top, etc. # One way is to define this permutation in 0-based indexing: # row0 --> row2, row1 --> row1, row2 --> row3, row3 --> row0 # or you can read it from your final P in the problem statement. # The user’s final P matrix is: # [0 0 0 1] # [0 1 0 0] # [1 0 0 0] # [0 0 1 0] # which sends row0->row2, row1->row1, row2->row3, row3->row0 in 0-based form. # Let’s encode that explicitly: perm = [2, 1, 3, 0] # meaning new row0=old row2, row1=old row1, etc. # Build P from this 'perm' list P = np.zeros((n, n), dtype=float) for new_row, old_row in enumerate(perm): P[new_row, old_row] = 1.0 # Now form the matrix Ap = P*A Ap = P @ A # Factor Ap by unpivoted (outer-product) LU U = Ap.copy() L = np.eye(n, dtype=float) for k in range(n - 1): pivot_val = U[k, k] # Fill in the multipliers below pivot_val for i in range(k+1, n): L[i, k] = U[i, k] / pivot_val # Rank-1 update on row i U[i, k+1:] = U[i, k+1:] - L[i, k] * U[k, k+1:] return L, U, P # -------------- Demo -------------- if __name__ == "__main__": A_test = np.array([ [ 4, 3, 2, 1], [ 8, 8, 5, 2], [16, 12, 10, 5], [32, 24, 20, 11] ]) L, U, P = fixedSolutionLU(A_test) print("L =\n", L, "\n") print("U =\n", U, "\n") print("P =\n", P, "\n") # Check that P*A == L*U: lhs = P @ A_test rhs = L @ U print("P*A =\n", lhs, "\n") print("L*U =\n", rhs, "\n") diff = np.linalg.norm(lhs - rhs) print(f"||P*A - L*U|| = {diff:e}")
Editor is loading...
Leave a Comment