Solving System of linear equations
Solving System of Linear Equations¶
- Gaussian Elimination
- LU Decomposition
- Inverse Matrix
- Gauss-Sidel Iteration
- Jecobi interation Method
Gaussian Elimination¶
In gaussian elimination method we find a row achilion matrix for the equation and just do back susbtitution.
import numpy as np
B = np.array([14.,13.,17.])
M = np.array([[1.,5.,1.],[2.,1.,3.],[3.,1.,4.]])
def gaussianElimination(M,B):
n = len(B)
x = np.zeros(n)
for i in range(n-1):
if M[i][i] == 0:
M[i+1], M[i] = M[i], M[i+1]
for j in range(i+1,n):
fector = M[j][i]/M[i][i]
for k in range(i,n):
M[j][k] += -M[i][k]*fector
B[j] -= B[i]*fector
#backsubstitution
for i in range(n-1,-1,-1):
s = 0
for j in range(i,n):
s += M[i][j]*x[j]
x[i] = (B[i]-s)/M[i][i]
return x
gaussianElimination(M,B)
LU Decomposition¶
Consdier a set of linear equations represented by the the matrix equation
$$MX = B$$Where $M$ is called the cofficient matrix and $B$ is called Constant Matrix.
We can decompose this matrix $M$ into two matrix $L$ and $L$ where L is a lower triangular matrix and U is an upper triangular matrix. Such that
$$ M = LU $$Now putting back into the euqation we get $$ LUX = B $$
We can see this as two sperate diffrential equations as $$LY = B$$ and $$UX = Y$$
Because $L$ and $U$ are triangular matrices, these eqution leads to just substiution. We first get $Y$ by forth susbtitution with $L$ and then we do back susbtitution with U and get X.
Getting L and U matrix:¶
The relation to get U and L matrix are given below
$$U_{ij} = a_{ij} - \sum_{k=0}^{i-1}L_{ik}U_{kj} $$$$ L_{ij} =\frac{1}{U_{jj}}( a_{ij} - \sum_{k=0}^{j-1}L_{ik}U_{kj}) $$If we look for this for a moment, it might seem complicated because L depends on U and U depand on L. The way out of this is suggested by Dr. Gajendra Purohit. We calculate first row of U matrix because for first row we don't need the sum. Then we calculate the first column of L matrix using it. The we calculate second row of U matrix, here we need the sum and we have already calculted the needed part. Then we calculate second column of L matirx. Hence one row of U and one column of L going on..
#a function to get L and U matrix seperately
def LU_decomposition(M,B):
#get number of unknowns
n = len(B)
#<---variables--->
L = np.zeros((n,n))
U = np.zeros((n,n))
for i in range(n):
#getting the U matrix
for j in range(i,n):
s = 0 #dummy var for sum
for k in range(i):
s += L[i][k]*U[k][j]
U[i][j] = M[i][j] -s
#getting L Matrix
for j in range(i,n):
s = 0
for k in range(i):
s += L[j][k]*U[k][i]
L[j][i] = (M[j][i]-s)/U[i][i]
return L,U
def LU_solver(M,B):
n = len(B)
#get L & U matrix
L,U = LU_decomposition(M,B)
#<--variables-->
y = np.zeros(n)
x = np.zeros(n)
#forth subsitution
for i in range(n):
s = 0
for j in range(i):
s += L[i][j]*y[j]
y[i] = B[i]-s
#back substitution
for i in range(n-1,-1,-1):
s = 0
for j in range(n-1,i,-1):
s += U[i][j]*x[j]
x[i] = (y[i]-s)/U[i][i]
return x
LU_solver(M,B)
Jecobi Iteration Method¶
For Diagonally Dominant Matrix, we use iteration methods.In words, a diagonally dominant matrix is a square matrix such that in each row, the absolute value of the term on the diagonal is greater than or equal to the sum of absolute values of the rest of the terms in that row. A strictly diagonally dominant matrix is non-singular, i.e. has an inverse.
for example:
$$ \begin{bmatrix} 3 &-2 &1 \\ 1&-3 & 2 \\ -1& 2 & 4 \end{bmatrix} $$M = [[2,-1,0],[-1,2,-1],[0,-1,2]]
B = [7,1,1]
Consider a set of linear equations represented in matrix form
\begin{equation} \begin{bmatrix} a_{00} & a_{01} & a_{02} \\ a_{10} & a_{11} & a_{12} \\ a_{20} & a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x2 \end{bmatrix} = \begin{bmatrix} b_0 \\ b_1 \\ b_2 \end{bmatrix} \end{equation}This can aslo be written as $$a_{00} x_0 + a_{01} x_1 + a_{02} x_2 = b_0 $$ $$a_{01} x_0 + a_{01} x_1 + a_{01} x_2 = b_1 $$ $$a_{02} x_0 + a_{02} x_1 + a_{02} x_2 = b_2 $$
More Explicitly:
$$x_0 = \frac{1}{a_{00}}(b_0 - a_{01} x_1 - a_{02}x_2)$$ $$x_1 = \frac{1}{a_{11}}(b_1 - a_{03} x_3 - a_{02}x_2)$$$$x_2 = \frac{1}{a_{22}}(b_2 - a_{01} x_1 - a_{03}x_3)$$In a nice form
$$ x_0 = \frac{1}{a_{00}}(b_0 - \sum_{j=0 , j != 0}^{3} a_{0j}x_{j})$$$$ x_1 = \frac{1}{a_{11}}(b_1 - \sum_{j=0 , j != 1}^{3} a_{1j}x_{j})$$$$ x_2 = \frac{1}{a_{22}}(b_2 - \sum_{j=0 , j != 2}^{3} a_{2j}x_{j})$$For any n-dimenssional matrix
$$ x_i = \frac{1}{a_{ii}}(b_i - \sum_{j=i , j != i}^{n} a_{ij}x_{j})$$The method for Jecobi Iteration goes like this:
We start with an intitial guess $x^0 = [x^0_0, x^0_1, x^0_2] $ put in above equation one by one i.e we put $ x^0_1, x^0_2$ in an equation for $x_0$ and find a new value we call this $x^1_0$ (first Iteration) and put $ x^0_0, x^0_2$ in an equation for $x_1$ and find a new value and so on. Now we have a new set of xs... this is called an Iteration. We repeat the process taking the new values and find a sequence which converges to the solutions.
#first input M and B
def jecobi(M,B):
n = len(B) #number of unknowns
max_iter = 50
count = 0
x = np.zeros(n) #initial guess
for _ in range(max_iter):
# print(x)
temp = x
#calculate next
for i in range(n):
s = 0
for j in range(n):
if j != i:
s += M[i][j]*temp[j]
x[i] =(B[i]+s)/M[i][i]
return x
jecobi(M,B)
Gauss Siedel Iteration Method¶
Gauss Siedel Iteration Method is exactly like Jecobi. We start with an intitial guess $x^0 = [x^0_0, x^0_1, x^0_2] $ put in above equation one by one i.e we put $ x^0_1, x^0_2$ in an equation for $x_0$ and find a new value we call this $x^1_0$ (first Iteration) and unlike Jecobi Iteration we pu the new value $ x^1_0, x^0_2$ in an equation for $x_1$ .
def gauss_sidel(M,B):
n = len(B) #number of unknowns
max_iter = 50
count = 0
x = np.zeros(n) #initial guess
for _ in range(max_iter):
# print(x)
#calculate next
for i in range(n):
s = 0
for j in range(n):
if j != i:
s += M[i][j]*x[j]
x[i] =(B[i]+s)/M[i][i]
return x
gauss_sidel(M,B)
Case Study: finding an Inverse Matrix for a squre matrix¶
Let A be a squre matrix than B is called its inverse matrix if $$AB = I$$ where I is identity matrix.
we know that an identity Matrix is like this:
$$ \begin{bmatrix} 1 &0 &0 & \cdots &0&0 \\ 0 &1 &0 & \cdots &0&0 \\ 0 &0 &1 & \cdots &0&0 \\ \vdots &\vdots &\vdots & &\vdots&\vdots \\ 0 &0 &0 & \cdots &1&0 \\ 0 &0 &0 & \cdots &0&1 \end{bmatrix} $$So we an write the problem as set of marix equations $$AB_i = I_i$$
Where $B_i$ is $i_{th}$ column matrix of B and $I_i$ is a column matrix such that all of elements are zero except $i^{th}$.
SO our problem is just to find all $B_i$s and glue them together.
def inverse(M):
n = len(M)
B = np.zeros((n,n))
for i in range(n):
I_i = np.identity(n)[i]
xs = LU_solver(M, I_i)
for j in range(n):
B[j][i] = xs[j]
return B
B = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
B = inverse(M)
print(M)
print(B)
def matmul(A,B):
rA = len(A) #number of rows in A
rB = len(B) #number of rows in B
cA = len(A[0]) #number of columns in A
cB = len(B[0]) #number of columns in B
AB = np.zeros((rA, cB)) #result matrix
if cA == rB:
for i in range(rA):
for j in range(cB):
s = 0
for ii in range(rA):
s += A[i][ii]*B[j][ii]
AB[i][j] = s
return AB
print(matmul(M, B))
Comments
Post a Comment