pythonPivot

import numpy as np

A0 = np.array([[3, 1, 4, -1],
               [2, -2, -1, 2],
               [5, 7, 14, 8],
               [1, 3, 2, 4]])

b0 = np.array([7, 1, 20, -4])

A = A0.copy()
b = b0.copy()
n = A.shape[0]
I = np.eye(n)
k = 0  # Initialize phase counter

print("============Phase 1===============")

k = k + 1  # Reset/update phase counter
# Locate position of element largest modulus in the pivot element
r = k + np.argmax(np.abs(A[k-1:, k-1]))
# Create pivot matrix to interchange rows r & k
P = np.eye(n)
P[[k-1, r], :] = P[[r, k-1], :]
# Interchange rows r & k of A & b
A = np.dot(P, A)
b = np.dot(P, b)
# Create the Gaussian Transformation matrix
i = np.arange(k, n)
M1 = np.eye(n)
M1[i, k-1] = -A[i, k-1] / A[k-1, k-1]
# Perform Elementary row reduction
A = np.dot(M1, A)
b = np.dot(M1, b)

print(M1)
print(A)
print(b)

print("=====================End of Phase 1===============")

print("============Phase 2===============")
k = k + 1  # Reset/update phase counter
# Locate position of element largest modulus in the pivot element
r = k + np.argmax(np.abs(A[k-1:, k-1]))
# Create pivot matrix to interchange rows r & k
P = np.eye(n)
P[[k-1, r], :] = P[[r, k-1], :]
# Interchange rows r & k of A & b
A = np.dot(P, A)
b = np.dot(P, b)
# Create the Gaussian Transformation matrix
i = np.arange(k, n)
M2 = np.eye(n)
M2[i, k-1] = -A[i, k-1] / A[k-1, k-1]
# Perform Elementary row reduction
A = np.dot(M2, A)
b = np.dot(M2, b)

print(M2)
print(A)
print(b)

print("=====================End of Phase 2===============")

print("============Phase 3===============")
k = k + 1  # Reset/update phase counter
# Locate position of element largest modulus in the pivot element
r = k + np.argmax(np.abs(A[k-1:, k-1]))
# Create pivot matrix to interchange rows r & k
P = np.eye(n)
P[[k-1, r], :] = P[[r, k-1], :]
# Interchange rows r & k of A & b
A = np.dot(P, A)
b = np.dot(P, b)
# Create the Gaussian Transformation matrix
i = np.arange(k, n)
M3 = np.eye(n)
M3[i, k-1] = -A[i, k-1] / A[k-1, k-1]
# Perform Elementary row reduction
A = np.dot(M3, A)
b = np.dot(M3, b)

print(M3)
print(A)
print(b)

print("=====================End of Phase 3===============")

print('Finally, the upper triangular matrix, U, and the right-hand side vector, b, are given as')
U = A  # Upper triangular matrix

print('Upper triangular matrix, U')
print(U)

print('Right-hand side, b')
print(b)
print("==============End of Alternative #1===========")

print("Alternative #2: Using A Script")

A = A0.copy()
b = b0.copy()
n = A.shape[0]
I = np.eye(n)
S = np.max(np.abs(A), axis=1)
k = 0  # Initialize phase counter

for k in range(1, n):
    # Locate position of element largest modulus in the pivot element
    r = k + np.argmax(np.abs(A[k-1:, k-1]) / S[k-1:])
    # Permutation matrix to swap rows r & k
    Pk = np.eye(n)
    Pk[[k-1, r], :] = Pk[[r, k-1], :]
    # Swap rows r & k of A & b
    A = np.dot(Pk, A)
    b = np.dot(Pk, b)
    S = np.dot(Pk, S)
    # Create kth Gaussian Transformation matrix
    i = np.arange(k, n)
    Mk = np.eye(n)
    Mk[i, k-1] = -A[i, k-1] / A[k-1, k-1]
    # Perform Elementary row reduction
    A = np.dot(Mk, A)
    b = np.dot(Mk, b)

print("Finally, the upper triangular matrix, U, and the right-hand side vector, b are given as")
U = A  # Upper triangular matrix

print('Upper triangular matrix, U')
print(U)

print('Right-hand side, b')
print(b)

print("Solution Of Upper Triangular System Using Back Substitution Method")

print("Alternative #1: Step-By-Step Approach")
x = np.zeros(n)  # Preallocate Solution

print("Solve for the last (4th) unknown")
i = n  # Initialize row to the last row
x[i-1] = b[i-1] / U[i-1, i-1]
print(x[i-1])

print("Solve for the last but one (3rd) unknown")
i = n - 1
k = i + 1
x[i-1] = (b[i-1] - np.dot(U[i-1, k-1:], x[k-1:])) / U[i-1, i-1]
print(x[i-1])

print("Solve for the last but two (2nd) unknown")
i = n - 2
k = i + 1
x[i-1] = (b[i-1] - np.dot(U[i-1, k-1:], x[k-1:])) / U[i-1, i-1]
print(x[i-1])

print("Solve for the 1st unknown")
i = n - 3
k = i + 1
x[i-1] = (b[i-1] - np.dot(U[i-1, k-1:], x[k-1:])) / U[i-1, i-1]

print('Solution')
print(x)

print("Alternative #1: Using A Script")
x = np.zeros(n)  # Preallocate Solution

x[-1] = b[-1] / U[-1, -1]
for i in range(n-2, -1, -1):
    k = i + 1
    x[i] = (b[i] - np.dot(U[i, k:], x[k:])) / U[i, i]

print('Solution')
print(x)

x0 = np.linalg.solve(A0, b0)