# Solving Linear Systems: Direct Methods
<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a><br />This notebook by Xiaozhou Li is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  
All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).


```python
%matplotlib inline
import numpy as np
from scipy.linalg import lu
```

## Gaussian Elimination

The simplest version of Gaussian Elimination only involves two operations as follows:
* Add or substrct a muliple of one equation from another.
* Multiply an equation by a nonzero constant.

### Solving a linear system by Gaussian elimination
The general form of a linear system for $n$ equations in $n$ unknown can be wriiten as
\begin{equation}
\begin{pmatrix}{a_{11}} & {a_{12}} & {\cdots} & {a_{1 n}} \\ {a_{21}} & {a_{22}} & {\cdots} & {a_{2 n}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {a_{n 1}} & {a_{n 2}} & {\cdots} & {a_{n n}}\end{pmatrix}\begin{pmatrix}{x_{1}} \\ {x_{2}} \\ {\vdots} \\ {x_{n}}\end{pmatrix}=\begin{pmatrix}{b_{1}} \\ {b_{2}} \\ {\vdots} \\ {b_{n}}\end{pmatrix}
\end{equation}
  
To solve this linear system by Gaussian elimination 
1. Using the allowed row operations to eliminate the system to an upper triangular system
\begin{equation}
\begin{pmatrix}{a_{11}} & {a_{12}} & {\cdots} & {a_{1 n}} & {b_{1}} \\ {} & {a_{22}^{(1)}} & {\cdots} & {a_{2 n}^{(1)}} & {b_{2}^{(1)}} \\ {} & {} & {\ddots} & {\vdots} & {\vdots} \\ {} & {} & {} & {a_{n n}^{(n-1)}} & {b_{n}^{(n-1)}}\end{pmatrix}
\end{equation}
2. Using the back substituion (backsolving) to solve the upper triangular system, 
\begin{align}
x_{n}=& \frac{b_{n}^{(n-1)}}{a_{n-1}^{(n-1)}}  \\ x_{k}=& \frac{b_{k}^{(k-1)}-\sum_{j=k+1}^{n} a_{k j}^{(k-1)} x_{j} }{a_{k k}^{(k-1)}}, \quad k=n-1, \ldots, 1. 
\end{align}
Here, assuming $a^{(k-1)}_{kk} \neq 0,\,\, k = 1,\ldots,n$

**Example**

Using Gaussian elimination to solve the following linear system
\begin{align}
 10^{-20}x_1 + x_2 & = 1 \\
 x_1 + 2x_2 & = 4
\end{align}
The argumented matrix reads
\begin{pmatrix}
10^{-20} & 1 & 1 \\
1 & 2 & 4
\end{pmatrix}

1. Elimination: $\text{Row}_2 = \text{Row}_2 - 10^{20}\times\text{Row}_1$ 


```python
A = np.array([[1e-20, 1, 1], [1, 2, 4]])
A[1] = A[1] - 10**20*A[0]
print (A[1])
```

    [ 0.e+00 -1.e+20 -1.e+20]


2. The echelon form of the argumented matrix reads
\begin{pmatrix}
10^{-20} & 1 & 1 \\
0 & -10^{20} & -10^{20}
\end{pmatrix}
3. Using the back substituion, the solution is 
$$ x_2 = 1,\qquad x_1 = 0.$$

## The $PA = LU$ factorization
### Partial pivoting
At the start of classical Gaussian elimination of $n$ equations in $n$ unknowns, the first step is to use the diagonal element $a_{11}$ as a pivot to eliminate the first column.  The partial pivoting protocol consists of comparing numbers before carrying out each elimination step. The largest entry of the first column is located, and its row is swapped with the pivot row, in this case the top row.

**Example**

Solve the following linear system
\begin{align}
 10^{-20}x_1 + x_2 & = 1 \\
 x_1 + 2x_2 & = 4
\end{align}

The argumented matrix reads
\begin{pmatrix}
10^{-20} & 1 & 1 \\
1 & 2 & 4
\end{pmatrix}

1. Interchange $\text{Row}_1$ with $\text{Row}_2$
\begin{pmatrix}
1 & 2 & 4 \\
10^{-20} & 1 & 1 
\end{pmatrix}
2. Elimination: $\text{Row}_2 = \text{Row}_2 - 10^{-20}\times\text{Row}_1$ 


```python
A = np.array([[1, 2, 4], [1e-20, 1, 1]])
A[1] = A[1] - 10**(-20)*A[0]
print (A[1])
```

    [0. 1. 1.]


3. The echelon form of the argumented matrix reads
\begin{pmatrix}
1 & 2 & 4 \\
0 & 1 & 1
\end{pmatrix}
4. Using the back substituion, the solution is 
$$ x_2 = 1,\qquad x_1 = 2.$$

* Solving this equation by the build-in linear solver in numpy:


```python
A = np.array([[1e-20, 1], [1, 2]])
b = np.array([1, 4])

x = np.linalg.solve(A, b)
print (x)
```

    [2. 1.]


### Implementation
#### The $LU$ factorization


```python
def LU_factor(A):
    A = A.astype(np.float)
    print (A)
    n = len(A)
    for j in range(0,n-1):
        if A[j,j] != 0:
            for i in range(j+1,n):
                lam = A[i,j]/A[j,j]
                A[i,j+1:n] = A[i,j+1:n] - lam*A[j,j+1:n]
                A[i,j] = lam
    return A
```

**Example**
Find the LU factorization of matrix
$$
A=\begin{pmatrix}{2} & {4} & {4} & {2} \\ {3} & {3} & {12} & {6} \\ {2} & {4} & {-1} & {2} \\ {4} & {2} & {1} & {1}\end{pmatrix}
$$
and 
$$
A=\begin{pmatrix}{1} & {1} & {3} & {4} \\ {2} & {4} & {1} & {3} \\ {3} & {2} & {1} & {1} \\ {2} & {1} & {2} & {1}\end{pmatrix}
$$


```python
A = np.array([[2, 4, 4, 2], [3, 3, 12, 6], [2, 4, -1, 2], [4, 2, 1, 1]])
print (LU_factor(A))
```

    [[ 2.  4.  4.  2.]
     [ 3.  3. 12.  6.]
     [ 2.  4. -1.  2.]
     [ 4.  2.  1.  1.]]
    [[ 2.   4.   4.   2. ]
     [ 1.5 -3.   6.   3. ]
     [ 1.  -0.  -5.   0. ]
     [ 2.   2.   3.8 -9. ]]



```python
A = np.array([[1, 1, 3, 4], [2, 4, 1, 3], [3, 2, 1, 1], [2, 1, 2, 1]])
print (LU_factor(A))
```

    [[1. 1. 3. 4.]
     [2. 4. 1. 3.]
     [3. 2. 1. 1.]
     [2. 1. 2. 1.]]
    [[  1.           1.           3.           4.        ]
     [  2.           2.          -5.          -5.        ]
     [  3.          -0.5        -10.5        -13.5       ]
     [  2.          -0.5          0.61904762  -1.14285714]]


#### The $PA = LU$ factorization

All of the techniques described so far are implemented in Python (Scipy). The most sophisticated form of Gaussian elimination we have discussed is the $PA=LU$ (or $A = PLU$) factorization. Scipy package has the command accepts a square coefficient matrix $A$ and returns $P$ , $L$, and $U$, such that $A = PLU$.


```python
from scipy.linalg import lu
A = np.array([[2, 1, 5], [4, 4, -4], [1, 3, 1]])
P, L, U = lu(A)
print (P)
print (L)
print (U)
```

    [[0. 0. 1.]
     [1. 0. 0.]
     [0. 1. 0.]]
    [[ 1.    0.    0.  ]
     [ 0.25  1.    0.  ]
     [ 0.5  -0.5   1.  ]]
    [[ 4.  4. -4.]
     [ 0.  2.  2.]
     [ 0.  0.  8.]]



```python
A = np.array([[2, 4, 4, 2], [3, 3, 12, 6], [2, 4, -1, 2], [4, 2, 1, 1]])
P, L, U = lu(A)
print (P)
print (L)
print (U)
```

    [[0. 0. 0. 1.]
     [0. 0. 1. 0.]
     [0. 1. 0. 0.]
     [1. 0. 0. 0.]]
    [[1.         0.         0.         0.        ]
     [0.5        1.         0.         0.        ]
     [0.75       0.5        1.         0.        ]
     [0.5        1.         0.41666667 1.        ]]
    [[ 4.     2.     1.     1.   ]
     [ 0.     3.    -1.5    1.5  ]
     [ 0.     0.    12.     4.5  ]
     [ 0.     0.     0.    -1.875]]



```python
print (np.dot(P, np.dot(L, U)) - A)
```

    [[0. 0. 0. 0.]
     [0. 0. 0. 0.]
     [0. 0. 0. 0.]
     [0. 0. 0. 0.]]


## Erros, Conditions Number

**Example** (Hilbert Matrix)

The $n\times n$ Hilbert Matrix $H_n$ is defined as follows:
$$
H_n = \begin{pmatrix}{1} & {1 / 2} & {1 / 3} & {\cdots} & {1 / n} \\ {1 / 2} & {1 / 3} & {1 / 4} & {\cdots} & {1 /(n+1)} \\ {1 / 3} & {1 / 4} & {1 / 5} & {\cdots} & {1 /(n+2)} \\ {\cdots} & {\cdots} & {\cdots} & {\cdots} & {\cdots} \\ {1 / n} & {1 /(n+1)} & {1 /(n+2)} & {\cdots} & {1 /(2 n-1)}\end{pmatrix}
$$

Solving the linear system  
$$ H_n x = b$$
with 
$$ b = H_n \cdot\begin{pmatrix}1 \\ 1 \\ \vdots \\ 1\end{pmatrix} $$
for $n = 5, 10, 20$
* The exact solution is $x = \begin{pmatrix}1 \\ 1 \\ \vdots \\ 1\end{pmatrix} $


```python
def hil(n):
    A = np.empty([n, n])
    for i in range(n):
        for j in range(n):
            A[i,j] = 1/(i+j+1)
    return A

def hil_example(n):
    A = hil(n)
    x_exact = np.zeros(n) + 1.
    b = np.dot(A, x_exact)
    x = np.linalg.solve(A, b)
    
    print ("The Hilbert Example")
    print ("Exact Solution: ", x_exact)
    print ("Numerical Solution: ", x)
    print ("Max Norm Error: ", np.max(x_exact - x))
```


```python
n = 30
#A = hil(n)
#print (A, '\n')

hil_example(n)
```

    The Hilbert Example
    Exact Solution:  [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
     1. 1. 1. 1. 1. 1.]
    Numerical Solution:  [  0.99999999   1.00000549   0.9995934    1.01063208   0.86492543
       1.95654753  -2.98980522  10.66985353 -10.64549422   1.40496892
      15.91216488  -9.67086839  -3.19672645   5.14709635  -5.92419006
      16.63359268   3.83136369 -13.36494379  -3.48351104  -2.37205148
      14.37850702  22.52283623 -37.62581185  20.90624367  -8.37575394
       1.62476071  10.16821532  -4.99742976   1.20141018   1.41386896]
    Max Norm Error:  38.62581184885538



```python
#n = 5
A = hil(n)
print (np.linalg.cond(A, np.inf))
```

    6.598338740631285e+18


**Example**

Solve the following linear system
\begin{align}
 (1 + 10^{-20})x_1 + x_2 & = 2 + 10^{-20} \\
 x_1 + x_2 & = 2
\end{align}


```python
A = np.array([[1e-20, 1], [1, 1]])
b = np.array([2+1e-20, 2])

x = np.linalg.solve(A, b)
print (x)
```

    [0. 2.]



```python

```

6e42003939665170c37eeff30f984ab972f5bbec
ipynb
Jupyter Notebook
LinearSys_DirectMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
2021-11-06T02:11:33.000Z
2021-11-06T02:11:33.000Z
LinearSys_DirectMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
LinearSys_DirectMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
Qwen/Qwen-72B
1. YES
2. YES