# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/dyjdlopez/numeth2021/blob/main/Week%207-9%20-%20Solving%20System%20of%20Linear%20Equations/NuMeth_3_System_of_Linear_Equations.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="pmwtuuOZ-40s"
# # System of Linear Equations
# $_{\text{©<NAME> | 2021 | Computational Methods for Computer Engineers}}$
#
# Another fundamental engineering computational technique is solving a system of linear equations. Linear systems are very much prevalent in engineering and computational modeling, and solving them is very crucial in analyzing and operating these systems. Solving systems of linear equations are useful in several practices such as optimization. Some non-engineering applications include finance, computer programming, geometry, etc. The coverage of the module is as follows:
# * Review of System of Linear Equations
# * Gaussian Elimintation
# * Gauss-Jordan Elimination
# * Cramer's Rule
# * Cholesky's LU Decomposition
# * Jacobi Method
# * Gauss-Siedel Method 
# * Successive over-relaxation (SOR) Method
# * Python Functions for Roots
# * Applications of Root-finding

# + [markdown] id="k_C7E7qJAy7K"
# ## 3.1 Looking back at Systems of Linear Equations
#
# A linear equation is described as an equation that has a scalar gradient. Or simply, it is an equation that produces a line when visualized. Formally,a  linear equation in $n$ (unknown) variables $x_1, ... , x_n$ has the form [[1]](https://mandal.ku.edu/math290/FTwelveM290/chpterOne.pdf):
# $$a_1x_1 + a_2x_2 + ... + a_nx_n = b 
# \\ _{\text{(Eq. 3.1.1)}}$$
# Whereas, $a_1, a_2,\dots,a_n,b$ are real where specifically $b$ is the constant term  and $a$ is the coefficient of $x_i$.
#
# And such, a system of these equations tell us that they are somewhat related or ar relevant for a specific scenario, event, or environment. Thus formally, $y$ a System of Linear Equations in n variables $x_1, x_2, ... , x_n$ we mean a collection of linear equations in these variables. A system of $m$ linear equations in these $n$ variables can be written as:
#
# $$\left\{
#     \begin{array}\\
#         a_{11}x_1 + a_{12}x_2 + ... + a_{1n}x_n = b_1\\ 
#         a_{21}x_1 + a_{22}x_2 + ... + a_{2n}x_n = b_2\\
#         a_{31}x_1 + a_{32}x_2 + ... + a_{3n}x_n = b_3\\
#         \vdots \\
#         a_{m1}x_1 + a_{m2}x_2 + ... + a_{mn}x_n = b_m\\
#     \end{array}
# \right.
# \\ _{\text{(Eq. 3.1.2)}}$$
# Whereas $a_{ij}$ and $b_i$ are all real numbers.
#
# If you can recall form Linear Algebra, Eq. 3.1.2 can be represented in augmented matrix form as: 
# $$\left[\begin{array}{cccc|c}
# a_{11}&a_{12}&\dots&a_{1n}&b_1\\
# a_{21}&a_{22}&\dots&a_{2n}&b_2\\
# a_{31}&a_{32}&\dots&a_{3n}&b_3\\
# \vdots&\vdots&\ddots&\vdots&\vdots\\
# a_{m1}&a_{m2}&\dots&a_{mn}&b_m
# \end{array}\right] \\ _{\text{(Eq. 3.1.3)}}$$
# Or in linear combination matrix form:
# $$\begin{bmatrix}
# a_{11}&a_{12}&\dots&a_{1n}\\
# a_{21}&a_{22}&\dots&a_{2n}\\
# a_{31}&a_{32}&\dots&a_{3n}\\
# \vdots&\vdots&\ddots&\vdots\\
# a_{m1}&a_{m2}&\dots&a_{mn}
# \end{bmatrix} \cdot
# \begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_n\end{bmatrix}=
# \begin{bmatrix}b_1\\b_2\\b_3\\\vdots\\b_m\end{bmatrix}
# \\ _{\text{(Eq. 3.1.4)}}$$

# + [markdown] id="Eb5N8gm0Jk4E"
# ### **Classifications of Linear Systems**
#
# Given a linear system in $n$ variables, precisely on the the following three is true [1]:
# * The system has **exactly one** solution (consistent system) — Geometrically, solution given by precisely the point where the graphs
# (two lines) of these two equations meet.
# * The system has **infinitely many** solutions (consistent system).
# * The system has **NO** solution (inconsistent system) — Geometrically, these two equations in the system represent two parallel lines (they never meet).
#
# Another classification exists for linear systems wherein:
# $$b_1=b_2=\dots=b_m=0 \\ _{\text{(Eq. 3.1.5)}}$$
# Which signifies the linear system is a **homogeneous linear system**.
#

# + [markdown] id="7Y3z7asyND8s"
# ### **Computational Methods**
# There are two types of computational methods in solving linear systems:
#
# * **Elimination Methods** - which eliminates parts of the coefficients matrix by means of elementary row and column operations in order to reach a form that yields the solution by simple calculations. Also known as direct method.
# * **Iterative Methods** - where equations are rearranged in a way that enables recursive calculations of values until convergence.

# + [markdown] id="VyH6BOT2K_aM"
# Now that you have a review on linear systems, we can now get started with solving them using computational methods.

# + id="jAnTIMT-JzS3"
import numpy as np

test_systems = {
    'test1': {
        'X': np.array([
                        [2, 7, -1, 3, 1],
                        [2, 3, 4, 1, 7],
                        [6, 2, -3, 2, -1],
                        [2, 1, 2, -1, 2],
                        [3, 4, 1, -2, 1]],float),
        'b': np.array([5, 7, 2, 3, 4], float),
        'ans': np.array([0.444, 0.556, 0.667, 0.222, 0.222], float)
    },
    'test2': {
        'X': np.array([
                        [0, 7, -1, 3, 1],
                        [2, 3, 4, 1, 7],
                        [6, 2, 0, 2, -1],
                        [2, 1, 2, 0, 2],
                        [3, 4, 1, -2, 1]],float),
        'b': np.array([5, 7, 2, 3, 4], float),
        'ans': np.array([0.021705, 0.792248, 1.051163, 0.158140, 0.031008], float)
},
'test3': {
        'X': np.array([
                        [2, -1, 5, 1],
                        [3, 2, 2, -6],
                        [1, 3, 3, -1],
                        [5, -2, -3, 3]],float),
        'b': np.array([5, 7, 2, 3, 4], float),
        'ans': np.array([2.0, -12.0, -4.0, 1.0], float)
},
'test4': {
        'X': np.array([
                        [8.00, 3.22, 0.80, 0.00, 4.10],
                        [3.22, 7.76, 2.33, 1.91, -1.03],
                        [0.80, 2.33, 5.25, 1.00, 3.02],
                        [0.00, 1.91, 1.00, 7.50, 1.03],
                        [4.10, -1.03, 3.02, 1.03, 6.44]],float),
        'b': np.array([9.45, -12.2, 7.78, -8.10, 10.0], float),
        'ans': np.array([2.0, -12.0, -4.0, 1.0], float)
}}

# + [markdown] id="Q564wD4cJ26q"
# ## 3.2 Gaussian Elimination
# The first method we will learn is the Gaussian Elimination. The goal of the Gaussian Elimination method is to reduce the linear system to its row-echelon form using basic or elementary row operations. Gaussian Elimination consists two stages: 
# * Elimination of the elements under the main diagonal of the coefficients matrix; and
# * Back-substitution of the solved unknowns until all system is solved.
#

# + id="iBtdBau1bIIR"
X = test_systems['test2']['X']
b = test_systems['test2']['b']
n = b.size
v = np.empty(n)

# + [markdown] id="t2HlC1L_Z0YG"
# ### Elimination
# 1. Loop of $k$ from $1$ to $n-1$ to index the fixed rows and eliminated columns
# 2. Loop of $i$ from $k+1$ to $n$ to index the subtracted rows
# 3. Loop of $j$ from $k$ to $n$ to index the columns for element subtraction the elimination statement using Eq. 3.2.1 and 3.2.2.
# $$x_{ij} = x_{kj} - \frac{x_{kk}}{x_{ik}}x_{ij} \\ _{\text{(Eq. 3.2.1)}}$$
# $$b_i = b_k - \frac{x_{kk}}{x_{ik}}b_i \\ _{\text{(Eq. 3.2.2)}}$$

# + id="9rJsipy8L9It"
for k in range(n-1):
    if X[k, k] == 0:
        for j in range (n):
            X[k,j], X[k+1, j] = X[k+1, j], X[k,j]
        b[k], b[k+1] = b[k+1], b[k]
    for i in range(k+1, n):
        if X[i, k] == 0: continue
        fctr = X[k, k] / X[i, k]
        b[i] = b[k] - fctr*b[i]
        for j in range(k, n):
            X[i, j] = X[k, j] - fctr*X[i, j]

# + [markdown] id="FUA2X7lraMWJ"
# ### Back-substitution
# 1. Starting with last row, let $v_n = \frac{b_n}{a_{nn}}$.
# 2. For rows from $i = n-1$ to $1$: substitute values of obtained for $v_i+1$ and compute the $v_i$ values of the current row. This can be expressed in the formula:
# $$v_i = \frac{b_i - \sum^n_{j=i+1}{a_{ij}v_j}}{a_{ii}} \\ _{\text{(Eq. 3.2.3)}}$$

# + id="J7P1kov-ZzFp" colab={"base_uri": "https://localhost:8080/"} outputId="769c30ae-1d97-4eaf-efb2-f9a8278d1766"
v[n-1] = b[n-1] / X[n-1, n-1]
for i in range(n-2, -1, -1):
    terms = 0
    for j in range(i+1, n):
        terms += X[i, j]*v[j]
    v[i] = (b[i] - terms)/X[i, i]
print('The solution of the system:')
print(v)


# + [markdown] id="JKCS8rbDSE8f"
# ### Activity 3.2
# Create a pseudocode and flowchart for the Gaussian Elimination method. You may use Eq. 3.2.1 to 3.2.3 for computational reference.

# + [markdown] id="sF-ZVORJL9k8"
# ## 3.3 Gauss-Jordan Elimination
# The Gauss-Jordan Elimination is quite similar to the Gauss elimination method. However, it differs in several aspects:
# * all elements above and below the main diagonal
# are eliminated to transform the coefficient matrix to the identity matrix $I_n$ and, accordingly, transform the constant terms vector to the solution vector.
# * for each pivot row, elimination operations are performed for all other rows above and below the pivot row. 
# * no back-substitution is applied since at the end of the elimination,
# the constant term vector will be already transformed to the solution vector.
#
# The following are the steps in performing Gauss-Jordan Elimination:
# 1. Performing partial pivoting by rearranging the rows of the system at each transformation to guarantee non-zero elements in the main diagonal of the coefficients matrix.
# 2. Division of every row of the system (including constant term) by the pivot element of that row ($x_{kk}$). This step will transform the matrix into an identity matrix.
# 3. Elimination of all elements above and below the main diagonal with applying the same operations on the constant terms vector. This algorithm requires the following three nested loops:
# >* The main loop of $k$ from row 1 to n to index the pivot rows. In this loop each element in the pivot row is divided by the pivot:
# $$x^*_{kj} = \frac{x_{kj}}{x_{kk}}, b^*_k = \frac{b_k}{x_{kk}}, k = 1,2,...,n, j=k,n$$
# >* The loop of $i$ from row $1$ to $n$ to index the subtraction rows. In order to avoid the subtraction of the pivot row from itself when $i$ is equal to $k$ or if the coefficient is already is equal to zero, the subtraction for current $i$ should be skipped.
# >*  The loop of j from k to n to index the columns for element subtraction starting from the pivot element to the last column of the coefficients matrix. Thus, the elimination statement is written as:
# $$x^*_{ij}=x_{ij}-x_{ik}x^*_{kj}$$
# and the new constant in the $i$ loop is:
# $$b^*_i=b_i-x_{ik}b^*_k$$
# Where, $k = 1,2,...,n i=1,2,...,n, i \not=k$.
#
#
#

# + id="8JEywa9IoWBg"
X = test_systems['test1']['X']
b = test_systems['test1']['b']
n = b.size
v = np.empty(n)


# + id="KZVh0JsdMApM"
def gssjrdn(X,b):
  a = np.array(X, float)
  b = np.array(b, float)
  n = len(b)
  for k in range(n):
    #Partial Pivoting
    if np.fabs(a[k,k]) < 1.0e-12:
      for i in range(k+1,n):
        if np.fabs(a[i,k]) > np.fabs(a[k,k]):
          a[[k,i]] = a[[i,k]]
          b[[k,i]] = b[[i,k]]
          break
    #Division of the pivot row
    pivot = a[k,k]
    a[k] /= pivot
    b[k] /= pivot
    #Elimination loop
    for i in range(n):
      if i == k or a[i,k] == 0: continue
      factor = a[i,k]
      a[i] -= factor * a[k]
      b[i] -= factor * b[k]
  return b,a


# + colab={"base_uri": "https://localhost:8080/"} id="0eX5WTDlG70l" outputId="bb34889f-c3a7-48e2-ebaf-1531dc41197c"
x,a = gssjrdn(X,b)

print("The solution of the system:")
print(x)
print("The coefficient matrix after transformation:")
print(a)

# + [markdown] id="o82lQAKqMAvM"
# ## 3.4 Cramer's Rule
# Is an explicit formula for the solution of a system of linear equations with as many equations as unknowns, valid whenever the system has a unique solution. It expresses the solution in terms of the determinants of the (square) coefficient matrix and of matrices obtained from it by replacing one column by the vector of right hand sides of the equations.
#
# 1. Set the coefficient matrix $X$ and the constants vector $b$.
# 2. Initialize an empty array of unknowns $v$ with a length of $m$.
# 3. For each element of $v$ their values are solved by:
# >* substituting the column of $X$ with $b$ corresponding to the index of the unknown element $v_i$, let this be matrix $C_i$.
# >* solve for $v_i$ by:
# $$v_i = \frac{|C_i|}{|X|}$$

# + id="81Etp-hQMJpV"
X = test_systems['test1']['X']
b = test_systems['test1']['b']


# + id="Zv9Ur4BOh4Ce"
def cramers(X,b):
  n = b.size
  v = np.empty(n)
  for i in range(b.size):
    C = X.copy()
    C[:,i] = b
    v[i] = np.linalg.det(C)/np.linalg.det(X)
  return v



# + colab={"base_uri": "https://localhost:8080/"} id="kL-lCZTlI422" outputId="c8e4a2a7-3244-4f02-c317-64a31d219a64"
cramers(X,b)

# + [markdown] id="DKGMJ3fWMJuU"
# ## 3.5 Cholesky Factorization
#
# LU decomposition of a matrix is the factorization of a given square matrix into two triangular matrices, one upper triangular matrix and one lower triangular matrix, such that the product of these two matrices gives the original matrix. It was introduced by Al<NAME> in 1948, who also created the turing machine [[2]](https://www.geeksforgeeks.org/l-u-decomposition-system-linear-equations/).
#
# Given a linear equation:
# $$AX=B\\_{\text{(Eq. 3.5.1)}}$$ it can be re-written as $$A = LU \\_{\text{(Eq. 3.5.2)}}$$ whereas $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. $L$ an $U$ are then called the decompositions of the matrix $A$. The decomposition can be illustrated as:
# $$\begin{bmatrix}
# a_{11}&a_{12}&a_{13}\\
# a_{21}&a_{22}&a_{23}\\
# a_{31}&a_{32}&a_{33}\\
# \end{bmatrix} =
# \begin{bmatrix}
# 1&0&0\\
# l_{21}&1&0\\
# l_{31}&l_{32}&1\\
# \end{bmatrix} \cdot
# \begin{bmatrix}
# u_{11}&u_{12}&u_{13}\\
# 0&u_{22}&u_{23}\\
# 0&0&u_{33}\\
# \end{bmatrix}
# \\ _{\text{(Eq. 3.5.3)}}$$
#
#
#
#
#

# + [markdown] id="UUSAWj4XTbAq"
# ### 3.5.1 Decomposition
# The **Cholesky decomposition** or **Cholesky factorization** is a decomposition of a symmetric and positive-definite matrix. The Cholesky decomposition is roughly twice as efficient as the LU decomposition for solving systems of linear equations.
#
# The Cholesky decomposition is a decomposition of the form $$A = LL^T$$ where $L$ is a lower triangular matrix with real and positive diagonal entries, and $L^T$ denotes the transpose of $L$.
#
# **Condition 1: Symmetry**
# To proceed with Cholesky's decomposition, the matrix should be symmetric. A matrix is said to be symmatric if it is a real square matrix $A$where:
# $$a_{ij}=a_{ji}, i=1,...,n, j=1,...,n\\ _{\text{(Eq. 3.5.4)}}$$ or simply, a matrix is symmetric if it is equal to its transpose. 
# $$A = A^T\\ _{\text{(Eq. 3.5.5)}}$$ 
#
# **Condition 2: Positive-Definite**
# A real symmetric matrix $A$ is positive definite if the eigenvalues of $A$ are positive.
#
# Once these conditions are met, we can then proceed with the Cholesky's Algorithm. The decomposed matrix can be solved by the following form:
# $$\begin{bmatrix}
# a_{11}& & & sym.\\
# a_{21}&a_{22}& & \\
# a_{31}&a_{32}&a_{33}& \\
# a_{41}&a_{42}&a_{43}&a_{44}\\
# \end{bmatrix} =
# \begin{bmatrix}
# l^2_{11}& & & sym.\\
# l_{11}l_{21}& l_{21}^2+l_{22}^2 & & \\
# l_{11}l_{31}& l_{31}l_{21} + l_{32}l_{22} &l_{31}^2+l_{32}^2+l_{33}^2& \\
# l_{11}l_{41}& l_{41}l_{21} + l_{42}l_{22} & l_{41}l_{31} + l_{42}l_{32} + l_{43}l_{33} &l_{41}^2+l_{42}^2+l_{43}^2+l_{44}^2 \\
# \end{bmatrix} 
# \\ _{\text{(Eq. 3.5.6)}}$$
#
# The algorithm can be shortened as:
# $$\left\{
#     \begin{array}\\
#         l_{ij} = \sqrt{a_{ij} - \sum_{k=1}^{j-1}l_{ik}^2} & i = j\\
#         l_{ij} = \frac{1}{l_{ij}} (a_{ij} - \sum_{k=1}^{j-1}l_{ik}l_{jk}) & i \not= j
#     \end{array}
# \right\}i=j,...,n; j=1,...,n 
# \\ _{\text{(Eq. 3.5.7)}}$$
#

# + id="OMJcJD65VmFM"
X = test_systems['test4']['X']
b = test_systems['test4']['b']


# + id="jGm4e_c-MUYk"
def cholesky(X):
  a = np.array(X, float)
  L = np.zeros_like(a)
  n = np.shape(a)[0]
  ## check for symmetry and positive-definite
  if (X == X.T).all() and (np.linalg.eigvals(X) > 0).all():
  
  ## Cholesky decomposition
    for j in range(n):
      for i in range(j,n):
        if i==j:
          L[i,j] = np.sqrt(a[i,j]-np.sum(L[i,:]**2))
        else:
          L[i,j]= (a[i,j]-np.sum(L[i,:j]*L[j,:j]))/L[j,j]
  else:
    print("Cannot perform Cholesky Factorization.")
  return L


# + colab={"base_uri": "https://localhost:8080/"} id="JketQe4JVvAm" outputId="b51480b4-a204-48f5-afe8-b32a6d1c2337"
cholesky(X)


# + [markdown] id="kZPcJxFnSxbL"
# ### 3.5.2 Forward Substitution
#
# $$\begin{bmatrix}
# l_{11}&0&0&... & 0\\
# l_{21}&l_{22}&0&...&0 \\
# l_{31}&l_{32}&l_{33}&...&0 \\
# \vdots&\vdots&\vdots&\ddots&\vdots\\
# l_{n1}&l_{n2}&l_{n3}&...&l_{nn}\\
# \end{bmatrix} 
# \left\{
#     \begin{array}\\
#         y_1\\y_2\\y_3\\\vdots\\y_n
#     \end{array}
# \right\}=
# \left\{
#     \begin{array}\\
#         b_1\\b_2\\b_3\\\vdots\\b_n
#     \end{array}
# \right\}
# \\ _{\text{(Eq. 3.5.8)}}$$
#
# The algorithm is represented as:
# $$y_i = \frac{1}{l_{ii}}(b_i - \sum^{i-1}_{j=1}l_{ij}y_j), i = 1,2,...,n\\ _{\text{(Eq. 3.5.9)}}$$

# + [markdown] id="vaS260i2VMzr"
# ### 3.5.3 Backward Substitution
# $$\begin{bmatrix}
# u_{1,1}&u_{1,2}&...&u_{1n-1}& u_{1,n}\\
# 0&u_{2,2}&...&u_{2,n-1}& u_{2,n} \\
# \vdots&\vdots&\vdots&\ddots&\vdots\\
# 0&0&...&u_{n-1,n-1}&u_{n-1,n}\\
# 0&0&...&0&u_{n,n}\\
# \end{bmatrix} 
# \left\{
#     \begin{array}\\
#         x_1\\x_2\\\vdots\\x_{n-1}\\x_n
#     \end{array}
# \right\}=
# \left\{
#     \begin{array}\\
#         y_1\\y_2\\\vdots\\y_{n-1}\\y_n
#     \end{array}
# \right\}
# \\ _{\text{(Eq. 3.5.10)}}$$
#
# The algorithm is represented as:
# $$y_i = \frac{1}{u_{ii}}(y_i - \sum^{n}_{j=i+1}u_{ij}x_j), i = n,n-1,...,1\\ _{\text{(Eq. 3.5.11)}}$$

# + id="I12DgFCzSoWM"
def solveLU(L,U,b):
  L = 

# + [markdown] id="G3cfZYs_MUSU"
# ## 3.6 Jacobi's Method

# + id="44J7r2y1MZy8"


# + [markdown] id="KmV2rF-zMZs8"
# ## 3.7 Gauss-Siedel Method

# + id="7q7x-KpFMeIU"


# + [markdown] id="TX_paBKtMdgd"
# ## 3.8 Successive over-relacation (SOR) Method 

# + id="L6y6Dlw4Mr3e"


# + [markdown] id="bUd1A-a_DVgi"
# ## References
# [1] <NAME>. (2020) *Introduction to Systems of Linear Equations*. Elementary Linear Algebra. Unversity of Kansas, KA, USA. https://mandal.ku.edu/math290/FTwelveM290/chpterOne.pdf
#

# + id="OfDLA99Nhixe"


# + id="Rdtbj_C_epg4"


# + id="BI_Tz6WGeq80"


StarcoderdataJupyterScriptsDedupFiltered