forked from tomer279/NumericalAnalysisMethods
-
Notifications
You must be signed in to change notification settings - Fork 0
/
GaussElimination.py
47 lines (43 loc) · 1.62 KB
/
GaussElimination.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
# Gaussian elimination with backward substitution
# to solve the linear equation Ax = b , such that A
# is a matrix of order n, and b is a column vector.
# The method receives matrix of order n x n+1, where the n+1 column
# represents the b vector
# the method returns a vector which is the solution for the linear equation.
# if there is no solution ,the method returns an appropriate print and stops.
# the method is based on algorithm 6.1 from "Numerical Analysis" by Burden, 10th Ed
# PARAMETERS:
# A -> a matrix of order n x (n+1), such that the first n columns create
# the coefficient matrix, and the last n+1 column is the vector b.
def gauss_elimination(A: np.matrix):
n = A.shape[0]
p = -1
for i in range(0, n - 1):
for p in range(i, n + 1):
if p == n:
print("no unique solution exists")
return
if A[p, i] != 0:
break
if p != i:
A[[p, i]] = A[[i, p]]
for j in range(i + 1, n):
A[j] = A[j] - (A[j, i] / A[i, i]) * A[i]
if A[n - 1, n - 1] == 0:
print("no unique solution exists")
return
x = np.zeros(n)
x[n - 1] = A[n - 1, n] / A[n - 1, n - 1]
for i in range(n - 1, -1, -1):
sig = 0
for j in range(i + 1, n):
sig += A[i, j] * x[j]
x[i] = (A[i, n] - sig) / A[i, i]
print(x)
return
A_mat = np.matrix([[1, 1, 0, 3, 4],
[2, 1, -1, 1, 1],
[3, -1, -1, 2, -3],
[-1, 2, 3, -1, 4]], dtype=float)
gauss_elimination(A_mat)