Using Python: In this problem, we\'ll do the LU decomposition. Again, you can ty
ID: 3808357 • Letter: U
Question
Using Python:
In this problem, we'll do the LU decomposition. Again, you can type in the system at the top of the file. Write a function that takes and n times n invertible matrix A and returns matrices P, L and U so that A = PLU where P is a permutation matrix, L is lower trianglar with 1's on the diagonal and U is upper trianbular. Use scaled partial pivoting. Write a function that takes a vector b and uses the decomposition PLU to solve the system Ax = b using backward substitution and forward substituion. Store the matrics as numpy arrays.Explanation / Answer
import pprint
import scipy
import scipy.linalg
A=scipy.array([7,3,-1,2],[3,8,1,-4],[-1,1,4,-1],[2,-4,-1,6])
P,L,U=scipy.linalg.lu(A)
Print ”A:”
Pprint.pprint(A)
Print “P:”
Pprint.pprint(P)
Print “L:”
Pprint.pprint(L)
Print “U:”
Pprint.pprint(U)
OUTPUT:-
Or
import numpy as np
from scipy.linalg import lu, inv
def gausselim(A,B):
"""
Solve Ax = B using Gaussian elimination and LU decomposition.
A = LU decompose A into lower and upper triangular matrices
LUx = B substitute into original equation for A
Let y = Ux and solve:
Ly = B --> y = (L^-1)B solve for y using "forward" substitution
Ux = y --> x = (U^-1)y solve for x using "backward" substitution
:param A: coefficients in Ax = B
:type A: numpy.ndarray of size (m, n)
:param B: dependent variable in Ax = B
:type B: numpy.ndarray of size (m, 1)
"""
# LU decomposition with pivot
pl, u = lu(A, permute_l=True)
# forward substitution to solve for Ly = B
y = np.zeros(B.size)
for m, b in enumerate(B.flatten()):
y[m] = b
# skip for loop if m == 0
if m:
for n in xrange(m):
y[m] -= y[n] * pl[m,n]
y[m] /= pl[m, m]
# backward substitution to solve for y = Ux
x = np.zeros(B.size)
lastidx = B.size - 1 # last index
for midx in xrange(B.size):
m = B.size - 1 - midx # backwards index
x[m] = y[m]
if midx:
for nidx in xrange(midx):
n = B.size - 1 - nidx
x[m] -= x[n] * u[m,n]
x[m] /= u[m, m]
return x
if __name__ == '__main__':
x = gausselim(np.array([[3, 2], [1, -4]]), np.array([[5], [10]]))
print x
OUTPUT:-
OR
import pprint
def mult_matrix(M, N):
"""Multiply square matrices of same dimension M and N"""
# Converts N into a list of tuples of columns
tuple_N = zip(*N)
# Nested list comprehension to calculate matrix multiplication
return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M]
def pivot_matrix(M):
"""Returns the pivoting matrix for M, used in Doolittle's method."""
m = len(M)
# Create an identity matrix, with floating point values
id_mat = [[float(i ==j) for i in xrange(m)] for j in xrange(m)]
# Rearrange the identity matrix such that the largest element of
# each column of M is placed on the diagonal of of M
for j in xrange(m):
row = max(xrange(j, m), key=lambda i: abs(M[i][j]))
if j != row:
# Swap the rows
id_mat[j], id_mat[row] = id_mat[row], id_mat[j]
return id_mat
def lu_decomposition(A):
"""Performs an LU Decomposition of A (which must be square)
into PA = LU. The function returns P, L and U."""
n = len(A)
# Create zero matrices for L and U
L = [[0.0] * n for i in xrange(n)]
U = [[0.0] * n for i in xrange(n)]
# Create the pivot matrix P and the multipled matrix PA
P = pivot_matrix(A)
PA = mult_matrix(P, A)
# Perform the LU Decomposition
for j in xrange(n):
# All diagonal entries of L are set to unity
L[j][j] = 1.0
# LaTeX: u_{ij} = a_{ij} - sum_{k=1}^{i-1} u_{kj} l_{ik}
for i in xrange(j+1):
s1 = sum(U[k][j] * L[i][k] for k in xrange(i))
U[i][j] = PA[i][j] - s1
# LaTeX: l_{ij} = rac{1}{u_{jj}} (a_{ij} - sum_{k=1}^{j-1} u_{kj} l_{ik} )
for i in xrange(j, n):
s2 = sum(U[k][j] * L[i][k] for k in xrange(j))
L[i][j] = (PA[i][j] - s2) / U[j][j]
return (P, L, U)
A = [[7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6]]
P, L, U = lu_decomposition(A)
print "A:"
pprint.pprint(A)
print "P:"
pprint.pprint(P)
print "L:"
pprint.pprint(L)
print "U:"
pprint.pprint(U)
OUTPUT:-
A:
[[7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6]]
P:
[[1.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 1.0]]
L:
[[1.0, 0.0, 0.0, 0.0],
[0.42857142857142855, 1.0, 0.0, 0.0],
[-0.14285714285714285, 0.2127659574468085, 1.0, 0.0],
[0.2857142857142857, -0.7234042553191489, 0.0898203592814371, 1.0]]
U:
[[7.0, 3.0, -1.0, 2.0],
[0.0, 6.714285714285714, 1.4285714285714286, -4.857142857142857],
[0.0, 0.0, 3.5531914893617023, 0.31914893617021267],
[0.0, 0.0, 0.0, 1.88622754491018]]
Related Questions
drjack9650@gmail.com
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.