numerik1/gaussian.py

76 lines
1.7 KiB
Python
Raw Normal View History

2018-11-11 17:49:43 +00:00
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)))