76 lines
1.7 KiB
Python
76 lines
1.7 KiB
Python
import copy
|
|
import pprint
|
|
import numpy
|
|
|
|
def gaussian_with_pivot(A, b, round = lambda x: x):
|
|
A = numpy.array(A, dtype = numpy.float64)
|
|
b = numpy.array(b, dtype = numpy.float64)
|
|
L = numpy.identity(len(A), dtype = numpy.float64)
|
|
s = numpy.zeros(len(A), dtype = numpy.float64)
|
|
|
|
for k in range(len(A)):
|
|
max_a = numpy.absolute(A[k:,k]).argmax()
|
|
A[[k, k + max_a]] = A[[k + max_a, k]]
|
|
b[[k, k + max_a]] = b[[k + max_a, k]]
|
|
s[k] = k + max_a
|
|
|
|
for i in range(k + 1, len(A)):
|
|
L[i][k] = round(A[i][k]/A[k][k])
|
|
b[i] = round(b[i] - round(L[i][k] * b[k]))
|
|
|
|
for j in range(k + 1, len(A)):
|
|
A[i][j] = round(A[i][j] - round(L[i][k] * A[k][j]))
|
|
|
|
A[i][k] = 0
|
|
return A, b, L, s
|
|
|
|
def gaussian_no_pivot(A, b, round = lambda x: x):
|
|
A = numpy.array(A, dtype = numpy.float64)
|
|
b = numpy.array(b, dtype = numpy.float64)
|
|
L = numpy.identity(len(A), dtype = numpy.float64)
|
|
|
|
for k in range(len(A)):
|
|
|
|
for i in range(k + 1, len(A)):
|
|
L[i][k] = round(A[i][k]/A[k][k])
|
|
b[i] = round(b[i] - round(L[i][k] * b[k]))
|
|
|
|
for j in range(k + 1, len(A)):
|
|
A[i][j] = round(A[i][j] - round(L[i][k] * A[k][j]))
|
|
|
|
A[i][k] = 0
|
|
return A, b, L
|
|
|
|
|
|
if( __name__ == "__main__"):
|
|
# For Exercise 1
|
|
A = [[10, 1, 0], [15, 2, 1], [5, 0, 1]]
|
|
b = [12, 18, 4]
|
|
pprint.pprint(gaussian_with_pivot(A, b))
|
|
|
|
|
|
# For Exercise 2
|
|
A = [[7, 1, 1]
|
|
, [10, 1, 1]
|
|
, [1000, 0, 1]]
|
|
b = [10, 13, 1001]
|
|
pprint.pprint(A)
|
|
|
|
def round4(x):
|
|
return float("%.4g" % x)
|
|
|
|
pprint.pprint(gaussian_with_pivot(A, b, round4))
|
|
pprint.pprint(gaussian_no_pivot(A, b, round4))
|
|
|
|
|
|
A = numpy.zeros((10, 10))
|
|
|
|
for i in range(10):
|
|
if(i - 1 > 0):
|
|
A[i][i - 1] = -1
|
|
A[i][i] = 2
|
|
if(i + 1 < 10):
|
|
A[i][i + 1] = -1
|
|
|
|
pprint.pprint(gaussian_no_pivot(A, numpy.ones(10)))
|