# Computational Astrophysics
## Linear Systems of Equations

---
## Eduard Larrañaga

Observatorio Astronómico Nacional\
Facultad de Ciencias\
Universidad Nacional de Colombia

---

### About this notebook

In this notebook we present some of the techniques to find the roots of a function.

---

A system of linear equations can be written in matrix form

\begin{equation}
A \mathbf{x} = \mathbf{b}
\end{equation}

where, $A$ is a real $n\times n$ matrix with coefficients $a_{ij}$.
$\mathbf{b}$ is a given real vector and $\mathbf{x}$ is a vector with
$n$ unknown values.


* If $\det{A} = |A| \ne 0$ and $\mathbf{b} \ne 0$, the  LSE has a unique solution
\begin{equation}
\mathbf{x} = A^{-1} \mathbf{b}\,\,,
\end{equation}
where $A^{-1}$ is the inverse of $A$ with $AA^{-1}=A^{-1}
A=\mathbf{I}$.


* If $\det{A} = 0$, the equations either
have no solution or an infinite number of solutions.



---
## Matrix Inversion Method

The inverse of a matrix $A$ is given by

\begin{equation}
A^{-1} = \frac{1}{|A|} \,\, \underbrace{\mathrm{adj}{A}}_\text{adjugate}\,\,.
\end{equation}

the adjugate of $A$ is the transpose of $A$'s cofactor matrix $C$:

\begin{equation}
\mathrm{adj}{A} = C^T\,\,.
\end{equation}

So the problem becomes finding $C$ and $\det A$. 


This method scales in complexity with $n!$, where n is the number of rows/columns of $A$.

--- 
## Upper-diagonal Matrix

Consider a linear system with an upper-diagonal matrix,

\begin{equation}
A\mathbf{x} = \begin{pmatrix}
a_{11}&a_{12}&a_{13}\\
0 &a_{22}&a_{23}\\
0 &0 &a_{33} \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}
= \begin{pmatrix}b_1\\b_2\\b_3\end{pmatrix}\,\,.
\end{equation}

This LSE is solved trivially by simple back-substitution:

\begin{equation}
x_3 = \frac{b_3}{a_{33}},\hspace{0.3cm}x_2 = \frac{1}{a_{22}}(b_2 - a_{23} x_3), \hspace{0.3cm} x_1 = \frac{1}{a_{11}} ( b_1 - a_{12} x_2 - a_{13} x_3).
\end{equation}


This is implemented in a simple way for a $n\times n$ system,


```python
import numpy as np

def backsubs(A,b):
    '''
    ------------------------------------------
    backsubs(A,b)
    ------------------------------------------
    Returns the solution fo the Linear System
    A x = b
    where
    A: nxn upper diagonal matrix
    b: n vector
    
    Arguments:
    A: upper-diagonal numpy array of size nxn
    b: numpy array of size n
    ------------------------------------------
    '''
    n = len(b)
    m = n-1
    # Create an empty vector for the solution
    x = np.zeros(n)
    
    # Main loop for back-substitution
    for i in range(m,-1,-1):
        xterms = 0
        for j in range(i+1,n):
            xterms = xterms + A[i,j]*x[j]
        x[i] = 1/A[i,i]*(b[i] - xterms)
    return x

```

**Example**

We will solve the LSE

\begin{equation}
\begin{pmatrix}
5 & 3 & 4 & 2\\
0 & -1 & 7 & 3 \\
0 & 0 & 5 & 4 \\
0 & 0 & 0 & 3 
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
x_4
\end{pmatrix} = 
\begin{pmatrix}
3 \\
14 \\
12 \\
6
\end{pmatrix}
\end{equation}


```python
A = np.array([[5., 3., 4., 2.], 
              [0., -1., 7., 3.], 
              [0., 0., 5., 4.], 
              [0., 0., 0., 3.]])

b = np.array([3., 14., 12., 6.])

x = backsubs(A,b)
x
```




    array([ 0.6, -2.4,  0.8,  2. ])



Lets probe that the solution is correct,


```python
np.dot(A,x)
```




    array([ 3., 14., 12.,  6.])



---
## Cramer's Rule

According to Cramer's rule the solution to the LSE is
 
\begin{equation}
x_i = \frac{\det A_i}{\det A}\,\,,
\end{equation}

where $A_i$ is the matrix formed from $A$ by replacing its $i$-th
column by the column vector $\mathbf{b}$.

Cramer's rule is more efficient than matrix inversion because it scales with $n^3$, so is more efficient
for large matrixes and has comparable efficiency to direct methods
such as Gauss Elimination.

### Example. 

Consider the system

\begin{equation}
\begin{pmatrix}
5&3&4\\
2&1&5\\
5&4&1
\end{pmatrix}
\begin{pmatrix} 
x_1\\ x_2\\ x_3  
\end{pmatrix}
= \begin{pmatrix}3\\4\\2\end{pmatrix}\,\,.
\end{equation}

The determinant of the matrix $A$ is

\begin{equation}
\det A = -14.
\end{equation}

Therefore we will use Cramer's Rule to solve. The solution for variable $x_1$ is given by 

\begin{equation}
x_1 = \frac{\det A_1}{\det A}\,\,,
\end{equation}

where $A_1$ is the matrix formed from $A$ by replacing its first
column by the column vector $\mathbf{b}$, i.e.

\begin{equation}
A_1 = \begin{pmatrix}
3&3&4\\
4&1&5\\
2&4&1
\end{pmatrix}.
\end{equation}

This matrix has $\det A_1 = 17$ and therefore

\begin{equation}
x_1 =- \frac{17}{14}.
\end{equation}

The solution for variable $x_2$ is given by 

\begin{equation}
x_2 = \frac{\det A_2}{\det A}\,\,,
\end{equation}

where $A_2$ is

\begin{equation}
A_2 = \begin{pmatrix}
5&3&4\\
2&4&5\\
5&2&1
\end{pmatrix}.
\end{equation}

This matrix has $\det A_2 = -25$ and therefore

\begin{equation}
x_2 = \frac{25}{14}.
\end{equation}

The solution for variable $x_3$ is given by 

\begin{equation}
x_3 = \frac{\det A_3}{\det A}\,\,,
\end{equation}

where $A_3$ is

\begin{equation}
A_3 = \begin{pmatrix}
5&3&3\\
2&1&4\\
5&4&2
\end{pmatrix}.
\end{equation}

This matrix has $\det A_3 = -13$ and therefore

\begin{equation}
x_3 = \frac{13}{14}.
\end{equation}

The complete solution is

\begin{equation}
\textbf{x}
= \frac{1}{14} \begin{pmatrix} -17 \\ 25 \\ 13 \end{pmatrix}\,\,.
\end{equation}

---

Now lets implement this example computationally.


```python
import numpy as np

# Defning matrix A
A = np.array([[5., 3., 4.], 
              [2., 1., 5.], 
              [5., 4., 1.]])

A
```




    array([[5., 3., 4.],
           [2., 1., 5.],
           [5., 4., 1.]])



We will use the function `numpy.linalg.det()` to calculate the determinants.


```python
detA = np.linalg.det(A)
detA
```




    -13.999999999999996



Since this determinant is non-zero, we define the vector $\textbf{b}$ and implement Cramer's rule,


```python
import CramersRule as cr
import numpy as np

A = np.array([[5., 3., 4.], 
              [2., 1., 5.], 
              [5., 4., 1.]])
b = np.array([3., 4., 2.])

x = cr.CramersRule(A,b)
x
```




    array([-1.21428571,  1.78571429,  0.92857143])



When the matrix is very large or when the values of the coefficients are very small, the function `numpy.linalg.det()` may give overflow or underflow errors. In these cases, it may be possible to calculate the determinant of the matrix using the function `numpy.linalg.slogdet()`. More information about this function at https://numpy.org/doc/stable/reference/generated/numpy.linalg.slogdet.html#numpy.linalg.slogdet

---
## Gauss Elimination

The Gauss algorithm consists of a series of steps to bring a $n \times n$ matrix into the  upper triangular form.

1. Sort the rows of $A$ so that the diagonal coefficient $a_{ii}$
  (called \emph{the pivot}) of row $i$ (for all $i$) is non-zero. If this
  is not possible, the LSE cannot be solved.


2. Replace the $j$-th equation with

\begin{equation}
-\frac{a_{j1}}{a_{11}} \times (\text{1-st equation}) + (j\text{-th equation})\,\,,
\end{equation}

where $j$ runs from $2$ to $n$. This will zero-out column 1 for $i>1$.

3. Repeat the previous step, but starting with the next row down and with 
$j >$ (current row number). The current row be row $k$. Then we must replace rows $j$, $j > k$, with 
  
  \begin{equation}
    -\frac{a_{jk}}{a_{kk}} \times (k\text{-th equation}) + (j\text{-th equation})\,\,,
  \end{equation}
  
where $k < j \le n$.

4. Repeat (3) until all rows have been reduced and the matrix is in upper
triangular form.


5. Back-substitute to find $\mathbf{x}$.


```python
import numpy as np
from GaussElimination import *

def backsubs(A,b):
    '''
    ------------------------------------------
    backsubs(A,b)
    ------------------------------------------
    Returns the solution fo the Linear System
    A x = b
    where
    A: nxn upper diagonal matrix
    b: n vector
    
    Arguments:
    A: upper-diagonal numpy array of size nxn
    b: numpy array of size n
    ------------------------------------------
    '''
    n = len(b)
    m = n-1
    # Create an empty vector for the solution
    x = np.zeros(n)
    
    # Main loop for back-substitution
    for i in range(m,-1,-1):
        xterms = 0
        for j in range(i+1,n):
            xterms = xterms + A[i,j]*x[j]
        x[i] = 1/A[i,i]*(b[i] - xterms)
    return x



A = np.array([[5., 3., 4.], [2., 1., 5.], [5., 4., 1.]])
b = np.array([3., 4., 2.])


A,b = GaussElim(A,b)

x = backsubs(A,b)
x
```




    array([-1.21428571,  1.78571429,  0.92857143])




```python
GaussElim?
```

---
## Tri-Diagonal Systems

Consider the following $4 \times 4 $ LSE, represented by a tri-diagonal matrix,

\begin{equation}
\begin{pmatrix}
b_1&c_1&0&0\\
a_1&b_2&c_2&0\\
0&a_2&b_3&c_3\\
0&0&a_3&b_4
\end{pmatrix}
\begin{pmatrix} 
x_1\\ x_2\\ x_3 \\ x_4 
\end{pmatrix}
= \begin{pmatrix}f_1\\f_2\\f_3\\f_4\end{pmatrix}\,\,.
\end{equation}


Tri-diagonal LSEs like this one are simple to solve using Gaussian elimination. When applying a forward elimination procedure, we obtain

\begin{equation}
\begin{pmatrix}
1&c_1/d_1&0&0\\
0&1&c_2/d_2&0\\
0&0&1&c_3/d_3\\
0&0&0&1
\end{pmatrix}
\begin{pmatrix} 
x_1\\ x_2\\ x_3 \\ x_4 
\end{pmatrix}
= \begin{pmatrix}y_1\\y_2\\y_3\\y_4\end{pmatrix}\,\,,
\end{equation}

where

\begin{equation}
\begin{array}{lcl}
d_1 &=& b_1\,,\\
d_2 &=& b_2 - a_{1} c_{1}/d_{1}\,,\\ 
d_3 &=& b_3 - a_{2} c_{2}/d_{2}\,,\\ 
d_4 &=& b_4 - a_{3} c_{3}/d_{3}\,,\\ 
\end{array}\hspace{2cm}
\begin{array}{lcl}
y_1 &=& f_1/d_1\,,\\
y_2 &=& (f_2 - y_1 a_1)/d_2\,,\\
y_3 &=& (f_3 - y_2 a_2)/d_3\,,\\
y_4 &=& (f_4 - y_3 a_3)/d_4\,.\\
\end{array}
\end{equation}

Then, using back-substitution, we obtain the solution

\begin{equation}
\begin{array}{lcl}
x_1 &=& y_1 - x_2 c_1 / d_1\,,\\
x_2 &=& y_2 - x_3 c_2 / d_2\,,\\
x_3 &=& y_3 - x_4 c_3 / d_3\,,\\
x_4 &=& y_4\,.
\end{array}
\end{equation}

--- 

In a general case, for a $n\times n$ system, we can write the solving algorithm as

1. **Forward Elimination.**
 * At the first step: $d_1 = b_1$ and $y_1 = f_1 / d_1$.
 
 * At the $k$-th step:
  \begin{equation*}
    \begin{aligned}
      d_k & = b_k - a_{k-1} c_{k-1} / d_{k-1}\,\,,\\
      y_k & = (f_k - y_{k-1} a_{k-1}) / d_k\,\,.
    \end{aligned}
  \end{equation*}

2. **Backward Substitution.** The $x_k$ are given by

  \begin{equation*}
    \begin{aligned}
      x_n &= y_n\,\,,\\
      x_{k-1} &= y_{k-1}- x_k c_{k-1} / d_{k-1}\,\,.
    \end{aligned}
  \end{equation*}



```python

```

7392fb2a41f67eaa096c59659e9038ea9ca97d20
ipynb
Jupyter Notebook
09._LSE/presentation/LSE01.ipynb
ashcat2005/ComputationalAstrophysics
edda507d0d0a433dfd674a2451d750cf6ad3f1b7
MIT
2020-09-23T02:49:10.000Z
2021-08-21T06:04:39.000Z
09._LSE/presentation/LSE01.ipynb
ashcat2005/ComputationalAstrophysics
edda507d0d0a433dfd674a2451d750cf6ad3f1b7
MIT
09._LSE/presentation/LSE01.ipynb
ashcat2005/ComputationalAstrophysics
edda507d0d0a433dfd674a2451d750cf6ad3f1b7
MIT
2020-12-05T14:06:28.000Z
2022-01-25T04:51:58.000Z
Qwen/Qwen-72B
1. YES
2. YES