Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy
- eps = 2**-52
- def LU(A):
- U = numpy.matrix(A)
- n = U.shape[0]
- assert len(U.shape)==2
- assert n==U.shape[1]
- L = numpy.matrix(numpy.identity(n))
- for j in range(n-1):
- assert abs(U[j,j])>=eps
- for i in range(j+1, n):
- mult=U[i,j]/U[j,j]
- for k in range(n):
- U[i,k] -= mult*U[j,k]
- L[i,j]=mult
- return L, U
- def backsolve(A, b):
- b = numpy.matrix(b).transpose()
- L, U = LU(A)
- n = L.shape[0]
- assert b.shape[0]==n
- assert b.shape[1]==1
- c = numpy.matrix(numpy.zeros((n, 1)))
- x = c.copy()
- for i in range(n):
- mult = c[i] = b[i]
- for k in range(i,n):
- b[k]-= mult*L[k,i]
- for i in range(n-1, -1, -1):
- mult = x[i] = c[i]/U[i,i]
- for k in range(i):
- c[k]-=mult*U[k,i]
- return x
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement