# CHAPTER 2 - Systems Of Equations


```python
# Import modules
import sys
import numpy as np
import numpy.linalg
import scipy
import sympy
import sympy.abc
from scipy import linalg
from scipy.sparse import linalg as slinalg
```

# 2.1 Gaussian Elimination


```python
def naive_gaussian_elimination(matrix):
    """
    A simple gaussian elimination to solve equations
    
    Args:
        matrix : numpy 2d array
        
    Returns:
        mat : The matrix processed by gaussian elimination
        x   : The roots of the equation
        
    Raises:
        ValueError:
            - matrix is null
        RuntimeError :
            - Zero pivot encountered
    """
    if matrix is None :
        raise ValueError('args matrix is null')
    
    #Clone the matrix
    mat = matrix.copy().astype(np.float64)
    
    # Row Size
    m = mat.shape[0]
    
    # Column Size
    n = mat.shape[1]
    
    # Gaussian Elimaination
    for i in range(0, m):
        if np.abs(mat[i , i]) == 0 :
            raise RuntimeError('zero pivot encountered')
        for j in range(i + 1, m):
            mult = mat[j, i] / mat[i, i]
            for k in range(i, m):
                mat[j, k] -= mult * mat[i, k]
            mat[j, n - 1] -= mult * mat[i, n - 1]
    
    # Back Substitution
    x = np.zeros(m, dtype=np.float64)
    for i in range(m - 1,-1,-1):
        for j in range(i + 1, m):
            mat[i, n-1] = mat[i ,n-1] - mat[i,j] * x[j]
            mat[i, j] = 0.0
        x[i] = mat[i, n-1] / mat[i, i]
        
    return mat, x
```

### Example

Apply Gaussian elimination in tableau form for the system of three equations in three
unknowns:

$$
\large
\begin{matrix}
x + 2y - z = 3 & \\ 
2x + y - 2z = 3 & \\ 
-3x + y + z = -6 
\end{matrix}
$$


```python
"""
Input:
[[ 1  2 -1  3]
 [ 2  1 -2  3]
 [-3  1  1 -6]]
"""
input_mat = np.array([1, 2, -1, 3, 2, 1, -2, 3, -3, 1, 1, -6])
input_mat = input_mat.reshape(3, 4)
output_mat, x = naive_gaussian_elimination(input_mat)

print(output_mat)
print('[x, y, z] = {}'.format(x))
```

    [[ 1.  0.  0.  3.]
     [ 0. -3.  0. -3.]
     [ 0.  0. -2. -4.]]
    [x, y, z] = [3. 1. 2.]


### Additional Examples

1. Put the system $x + 2y - z = 3,-3x + y + z = -6,2x + z = 8$ into tableau form and solve by Gaussian elimination.


```python
input_mat = np.array([
    [ 1, 2, -1,  3],
    [-3, 1,  1, -6],
    [ 2, 0,  1,  8]
])

output_mat, x = naive_gaussian_elimination(input_mat)

print(output_mat)
print('[x, y, z] = {}'.format(x))
```

    [[1.         0.         0.         3.        ]
     [0.         7.         0.         7.        ]
     [0.         0.         1.85714286 3.71428571]]
    [x, y, z] = [3. 1. 2.]


## 2.1 Computer Problems

1. Put together the code fragments in this section to create a MATLAB program for “naive” Gaussian elimination (meaning no row exchanges allowed). Use it to solve the systems of Exercise 2.

See my implementation **naive_gaussian_elimination** in python.

2. Let $H$ denote the $n \times n$ Hilbert matrix, whose $(i, j)$ entry is $1 / (i + j - 1)$. Use the MATLAB program from Computer Problem 1 to solve $Hx = b$, where $b$ is the vector of all ones, for (a) n = 2 (b) n = 5 (c) n = 10.


```python
def computer_problems2__2_1(n):
    # generate Hilbert matrix H
    H = scipy.linalg.hilbert(n)
    
    # generate b
    b = np.ones(n).reshape(n, 1)
    
    # combine H:b in tableau form
    mat = np.hstack((H, b))
    
    # gaussian elimination
    _, x = naive_gaussian_elimination(mat)
    
    return x
    
with np.printoptions(precision = 6, suppress = True):
    print('(a) n =  2 → x = {}'.format(computer_problems2__2_1( 2)))
    print('(b) n =  5 → x = {}'.format(computer_problems2__2_1( 5)))
    print('(c) n = 10 → x = {}'.format(computer_problems2__2_1(10)))
```

    (a) n =  2 → x = [-2.  6.]
    (b) n =  5 → x = [    5.  -120.   630. -1120.   630.]
    (c) n = 10 → x = [      -9.997365      989.771861   -23755.13378    240195.71429
     -1261048.597184  3783198.501116 -6725765.489567  7000357.237863
     -3937735.417591   923673.408496]


---
# 2.2 The LU Factorization


```python
def LU_factorization(matrix):
    """
    LU decomposition
    
    Arguments:
        matrix : numpy 2d array
        
    Return:
        L : lower triangular matrix
        U : upper triangular matrix
        
    Raises:
        ValueError:
            - matrix is null
            - matrix is not a 2d array
        RuntimeError :
            - zero pivot encountered
    """
    if matrix is None :
        raise ValueError('args matrix is null')
    
    if matrix.ndim != 2 :
        raise ValueError('matrix is not a 2d-array')
    
    # dimension
    dim = matrix.shape[0]
    
    # Prepare LU matrixs
    L = np.identity(dim).astype(np.float64)
    U = matrix.copy().astype(np.float64)
        
    # Gaussian Elimaination
    for i in range(0, dim - 1):
        # Check pivot is not zero
        if np.abs(U[i , i]) == 0 :
            raise RuntimeError('zero pivot encountered')
        for j in range(i + 1, dim):
            mult = U[j, i] / U[i, i]
            for k in range(i, dim):
                U[j, k] -= mult * U[i, k]
            L[j, i] = mult
        
    return L, U
```

## DEFINITION 2.2

An $m \times n$ matrix $L$ is **lower triangular** if its entries satisfy $l_{ij} = 0$ for $i < j$. An $m \times n$ matrix $U$ is **upper triangular** if its entries satisfy $u_{ij} = 0$ for $i > j$.

### Example

Find the LU factorization for the matrix $A$ in

$$
\large
\begin{bmatrix}
1 &  1 \\ 
3 & -4 \\ 
\end{bmatrix}
$$


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

L, U = LU_factorization(A)

print('L = ')
print(L)

print()

print('U = ')
print(U)
```

    L = 
    [[1. 0.]
     [3. 1.]]
    
    U = 
    [[ 1.  1.]
     [ 0. -7.]]


### Example

Find the LU factorization of A =

$$
\large
\begin{bmatrix}
1  & 2 & -1 \\ 
2  & 1 & -2 \\ 
-3 & 1 &  1 \\
\end{bmatrix}
$$


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

L, U = LU_factorization(A)

print('L = ')
print(L)

print()

print('U = ')
print(U)
```

    L = 
    [[ 1.          0.          0.        ]
     [ 2.          1.          0.        ]
     [-3.         -2.33333333  1.        ]]
    
    U = 
    [[ 1.  2. -1.]
     [ 0. -3.  0.]
     [ 0.  0. -2.]]


### Example

Solve system

$$
\large
\begin{bmatrix}
1 &  1 \\ 
3 & -4 \\ 
\end{bmatrix}
\begin{bmatrix}
x_1 \\ 
x_2 \\ 
\end{bmatrix}
=
\begin{bmatrix}
3 \\ 
2 \\ 
\end{bmatrix}
$$

, using the LU factorization


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

b = np.array([3, 2]).reshape(2, 1)

L, U = LU_factorization(A)

# calculate Lc = b where Ux = c
mat = np.hstack((L, b))
c = naive_gaussian_elimination(mat)[1].reshape(2, 1)

# calculate Ux = c
mat = np.hstack((U, c))
x = naive_gaussian_elimination(mat)[1].reshape(2, 1)

# output the result
print('x1 = {}, x2 = {}'.format(x[0][0], x[1][0]))
```

    x1 = 2.0, x2 = 1.0



### Example
Solve system
\begin{matrix}
x + 2y - z = 3 & \\ 
2x + y - 2z = 3 & \\ 
-3x + y + z = -6 
\end{matrix}
 
 using the LU factorization


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

b = np.array([3, 3, -6]).reshape(3, 1)

L, U = LU_factorization(A)

# calculate Lc = b where Ux = c
mat = np.hstack((L, b))
c = naive_gaussian_elimination(mat)[1].reshape(3, 1)

# calculate Ux = c
mat = np.hstack((U, c))
x = naive_gaussian_elimination(mat)[1].reshape(3, 1)

# output the result
print('x1 = {}, x2 = {}, x3 = {}'.format(x[0][0], x[1][0], x[2][0]))
```

    x1 = 3.0, x2 = 1.0, x3 = 2.0


### Additional Examples

1. Solve

$$
\large
\begin{bmatrix}
2 &  4 & -2 \\ 
1 & -2 &  1 \\ 
4 & -4 &  8 \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\ 
x_2 \\ 
x_3 \\
\end{bmatrix}
=
\begin{bmatrix}
6 \\ 
3 \\ 
0 \\
\end{bmatrix}
$$

using the A = LU factorization


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

b = np.array([6, 3, 0]).reshape(3, 1)

L, U = LU_factorization(A)

# calculate Lc = b where Ux = c
mat = np.hstack((L, b))
c = naive_gaussian_elimination(mat)[1].reshape(3, 1)

# calculate Ux = c
mat = np.hstack((U, c))
x = naive_gaussian_elimination(mat)[1].reshape(3, 1)

# output the result
print('x1 = {}, x2 = {}, x3 = {}'.format(x[0][0], x[1][0], x[2][0]))
```

    x1 = 3.0, x2 = -1.0, x3 = -2.0


## 2.2 Computer Problems

1. Use the code fragments for Gaussian elimination in the previous section to write a MATLAB script to take a matrix A as input and output L and U. No row exchanges are allowed - the program should be designed to shut down if it encounters a zero pivot. Check your program by factoring the matrices in Exercise 2.

See my implementation **LU_factorization** in python.


```python
# Exercise 2 - (a)
A = np.array([
    [ 3, 1, 2],
    [ 6, 3, 4],
    [ 3, 1, 5]
])

L, U = LU_factorization(A)

print('L = ')
print(L)

print()

print('U = ')
print(U)
```

    L = 
    [[1. 0. 0.]
     [2. 1. 0.]
     [1. 0. 1.]]
    
    U = 
    [[3. 1. 2.]
     [0. 1. 0.]
     [0. 0. 3.]]



```python
# Exercise 2 - (b)
A = np.array([
    [ 4, 2, 0],
    [ 4, 4, 2],
    [ 2, 2, 3]
])

L, U = LU_factorization(A)

print('L = ')
print(L)

print()

print('U = ')
print(U)
```

    L = 
    [[1.  0.  0. ]
     [1.  1.  0. ]
     [0.5 0.5 1. ]]
    
    U = 
    [[4. 2. 0.]
     [0. 2. 2.]
     [0. 0. 2.]]



```python
# Exercise 2 - (c)
A = np.array([
    [ 1, -1, 1,  2],
    [ 0,  2, 1,  0],
    [ 1,  3, 4,  4],
    [ 0,  2, 1, -1]
])

L, U = LU_factorization(A)

print('L = ')
print(L)

print()

print('U = ')
print(U)
```

    L = 
    [[1. 0. 0. 0.]
     [0. 1. 0. 0.]
     [1. 2. 1. 0.]
     [0. 1. 0. 1.]]
    
    U = 
    [[ 1. -1.  1.  2.]
     [ 0.  2.  1.  0.]
     [ 0.  0.  1.  2.]
     [ 0.  0.  0. -1.]]


2. Add two-step back substitution to your script from Computer Problem 1, and use it to solve the systems in Exercise 4.


```python
def LU_factorization_with_back_substitution(A, b):
    """
    LU decomposition with two-step back substitution
    where Ax = b
    
    Arguments:
        A : coefficient matrix
        b : constant vector
        
    Return:
        x : solution vector
    """
    L, U = LU_factorization(A)

    # row size
    rowsz = b.size
    
    # calculate Lc = b where Ux = c
    matrix = np.hstack((L, b))
    c = naive_gaussian_elimination(matrix)[1].reshape(rowsz, 1)

    # calculate Ux = c
    matrix = np.hstack((U, c))
    x = naive_gaussian_elimination(matrix)[1].reshape(rowsz)
        
    return x
```


```python
# Exercise 4 - (a)
A = np.array([
    [ 3, 1, 2],
    [ 6, 3, 4],
    [ 3, 1, 5]
])

b = np.array([0, 1, 3]).reshape(3, 1)

x = LU_factorization_with_back_substitution(A, b)

print(x)
```

    [-1.  1.  1.]



```python
# Exercise 4 - (b)
A = np.array([
    [ 4, 2, 0],
    [ 4, 4, 2],
    [ 2, 2, 3]
])

b = np.array([2, 4, 6]).reshape(3, 1)

x = LU_factorization_with_back_substitution(A, b)

print(x)
```

    [ 1. -1.  2.]


---
# 2.3 Sources Of Error

## DEFINITION 2.3

The **inﬁnity norm**, or **maximum norm**, of the vector $x = (x_1, \cdots, x_n)$ is $||x||_{\infty} = \text{max}|x_i|, i = 1,\cdots,n$, that is, the maximum of the absolute values of the components of x.

## DEFINITION 2.4

Let $x_a$ be an approximate solution of the linear system $Ax = b$. The **residual** is the vector $r = b - Ax_a$. The **backward error** is the norm of the residual $||b - Ax_a||_{\infty}$,
and the **forward error** is $||x - x_a||_{\infty}$.

### Example 

Find the backward and forward errors for the approximate solution $x_a = [1, 1]$ of the system

$$
\large
\begin{bmatrix}
1 &  1 \\ 
3 & -4 \\ 
\end{bmatrix}
\begin{bmatrix}
x_1 \\ 
x_2 \\ 
\end{bmatrix}
=
\begin{bmatrix}
3 \\ 
2 \\ 
\end{bmatrix}
$$


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

b = np.array([3, 2])

xa = np.array([1, 1])
```


```python
# Get correct solution
system = sympy.Matrix(((1, 1, 3), (3, -4, 2)))
solver = sympy.solve_linear_system(system, sympy.abc.x, sympy.abc.y)

# Packed as list
x = np.array([solver[sympy.abc.x].evalf(), solver[sympy.abc.y].evalf()])

# Output
print(x)
```

    [2.00000000000000 1.00000000000000]



```python
# Get backward error (differences in the input)
residual = b - np.matmul(A, xa)
backward_error = np.max(np.abs(residual))
print('backward error is {:f}'.format(backward_error))

# Get fowrawd error (differences in the output)
forward_error = np.max(np.abs(x - xa))
print('forward error is {:f}'.format(forward_error))
```

    backward error is 3.000000
    forward error is 1.00000000000000


### Example

Find the forward and backward errors for the approximate solution [-1, 3.0001] of the system

$$
\large
\begin{align*}
x_1 + x_2 &= 2 \\ 
1.0001 x_1 + x_2 &= 2.0001 \\
\end{align*}
$$


```python
A = np.array([
    [      1, 1],
    [ 1.0001, 1],
])

b = np.array([2, 2.0001])

# approximated solution
xa = np.array([-1, 3.0001])

# correct solution
x = LU_factorization_with_back_substitution(A, b.reshape(2, 1))
```


```python
# Get backward error 
residual = b - np.matmul(A, xa)
backward_error = np.max(np.abs(residual))
print('backward error is {:f}'.format(backward_error))

# Get fowrawd error 
forward_error = np.max(np.abs(x - xa))
print('forward error is {:f}'.format(forward_error))
```

    backward error is 0.000100
    forward error is 2.000100


The **relative backward error** of system $Ax = b$ is deﬁned to be $\large \frac{||r||_{\infty}}{||b||_{\infty}}$.

The **relative forward error** is $\large \frac{||x - x_a||_{\infty}}{||x||_{\infty}}$.

The **error magniﬁcation factor** for $Ax = b$ is the ratio of the two, or $\large \text{error magnification factor} = \frac{\text{relative forward error}}{\text{relative backward error}} = \frac{\frac{||x - x_a||_{\infty}}{||x||_{\infty}}}{\frac{||r||_{\infty}}{||b||_{\infty}}}$

## DEFINITION 2.5

The **condition number** of a square matrix A, **cond(A)**, is the maximum possible error magniﬁcation factor for solving Ax = b, over all right-hand sides b.

The matrix norm of an n x n matrix A as 

$$
\large ||A||_{\infty} = \text{maximum absolute row sum}
$$


```python
def matrix_norm(A):
    rowsum = np.sum(np.abs(A), axis = 1)
    return np.max(rowsum)
```

## THEOREM 2.6

The condition number of the n x n matrix A is

$$
\large cond(A) = ||A|| \cdot ||A^{-1}||
$$


```python
def condition_number(A):
    inv_A = np.linalg.inv(A)
    cond = matrix_norm(A) * matrix_norm(inv_A)
    return cond
```

### Additional Examples

1. Find the determinant and the condition number (in the inﬁnity norm) of the matrix

$$
\large
\begin{bmatrix}
811802 & 810901 \\ 
810901 & 810001 \\ 
\end{bmatrix}
$$


```python
A = np.array([
    [ 811802, 810901],
    [ 810901, 810001],
])

print('determinant of A is {}'.format(scipy.linalg.det(A)))
print('condition number : {:.4e}'.format(condition_number(A)))
```

    determinant of A is 1.0000645141117275
    condition number : 2.6330e+12


2. The solution of the system

$$
\large
\begin{bmatrix}
2 &  4.01 \\ 
3 & 6 \\ 
\end{bmatrix}
\begin{bmatrix}
x_1 \\ 
x_2 \\ 
\end{bmatrix}
=
\begin{bmatrix}
6.01 \\ 
9 \\ 
\end{bmatrix}
$$

is $[1, 1]$

(a) Find the relative forward and backward errors and error magniﬁcation (in the inﬁnity norm) for the approximate solution [21,-9].


```python
A = np.array([
    [ 2, 4.01],
    [ 3, 6.00],
])

b = np.array([6.01, 9])

# approximated solution
xa = np.array([21, -9])

# correct solution
x = LU_factorization_with_back_substitution(A, b.reshape(2, 1))

# forward error
forward_error = np.max(np.abs(x - xa))

# relative forward error
relative_forward_error = forward_error / np.max(np.abs(x))

# backward error
backward_error = np.max(np.abs(b - np.matmul(A, xa)))

# relative backward error
relative_backward_error = backward_error / np.max(np.abs(b))

# error magnification factor
error_magnification_factor = relative_forward_error / relative_backward_error

print('relative forward error : {}'.format(relative_forward_error))
print('relative backward error : {}'.format(relative_backward_error))
print('error magnification factor : {}'.format(error_magnification_factor))
```

    relative forward error : 19.999999999998934
    relative backward error : 0.011111111111110676
    error magnification factor : 1799.9999999999745


(b) Find the condition number of the coefﬁcient matrix.


```python
A = np.array([
    [ 2, 4.01],
    [ 3, 6.00],
])

print('condition number : {}'.format(condition_number(A)))
```

    condition number : 3003.000000000063


## 2.3 Computer Problems

1. For the n x n matrix with entries $A_{ij} = 5 / (i + 2j - 1)$, set $x = [1,\cdots,1]^T$ and $b = Ax$. Use the MATLAB program from Computer Problem 2.1.1 or MATLAB’s backslash command to compute $x_c$, the double precision computed solution. Find the inﬁnity norm of the forward error and the error magniﬁcation factor of the problem $Ax = b$, and compare it with the condition number of A: (a) n = 6 (b) n = 10.


```python
def system_provider(n, data_generator):
    A = np.zeros([n, n])
    x = np.ones(n)
    
    for i in range(n):
        for j in range(n):
            A[i, j] = data_generator(i + 1, j + 1)
    
    b = np.matmul(A, x)
            
    return A, x, b

def problem_2_3_1_generic_solver(n, data_generator):
    A, x, b = system_provider(n, data_generator)
    xc = np.linalg.solve(A, b)
    
    # forward error
    forward_error = np.max(np.abs(x - xc))
    
    # relative forward error
    relative_forward_error = forward_error / np.max(np.abs(x))
    
    # backward error
    backward_error = np.max(np.abs(b - np.matmul(A, xc)))
    
    # relative backward error
    relative_backward_error = backward_error / np.max(np.abs(b))
    
    # error magnification factor
    error_magnification_factor = relative_forward_error / relative_backward_error
    
    # condition number
    condA = condition_number(A)
    
    return forward_error, error_magnification_factor, condA
```


```python
def problem_2_3_1_solver(n):
    return problem_2_3_1_generic_solver(n, lambda i, j : 5 / (i + 2 * j - 1))

# (a) n = 6
print('(a) n =  6, forward error = {:.3g}, error magnification factor = {:.3g}, condition number = {:.3g}'.format(*problem_2_3_1_solver(6)))

# (b) n = 10
print('(b) n = 10, forward error = {:.3g}, error magnification factor = {:.3g}, condition number = {:.3g}'.format(*problem_2_3_1_solver(10)))
```

    (a) n =  6, forward error = 8.13e-10, error magnification factor = 5.61e+06, condition number = 7.03e+07
    (b) n = 10, forward error = 0.0017, error magnification factor = 1.4e+13, condition number = 1.31e+14


2. Carry out Computer Problem 1 for the matrix with entries $A_{ij} = 1/(|i - j| + 1)$.


```python
def problem_2_3_2_solver(n):
    return problem_2_3_1_generic_solver(n, lambda i, j : 1 / (np.abs(i - j) + 1))

# (a) n = 6
print('(a) n =  6, forward error = {:.3g}, error magnification factor = {:.3g}, condition number = {:.3g}'.format(*problem_2_3_2_solver(6)))

# (b) n = 10
print('(b) n = 10, forward error = {:.3g}, error magnification factor = {:.3g}, condition number = {:.3g}'.format(*problem_2_3_2_solver(10)))
```

    (a) n =  6, forward error = 4.44e-16, error magnification factor = 2.92, condition number = 8.61
    (b) n = 10, forward error = 1.22e-15, error magnification factor = 10.3, condition number = 11.3


3. Let A be the n x n matrix with entries $A_{ij} = |i - j| + 1$. Deﬁne $x = [1,\cdots,1]^T$ and $b = Ax$. For n = 100,200,300,400, and 500, use the MATLAB program from Computer Problem 2.1.1 or MATLAB’s backslash command to compute $x_c$, the double precision computed solution. Calculate the inﬁnity norm of the forward error for each solution. Find the ﬁve error magniﬁcation factors of the problems $Ax = b$, and compare with the corresponding condition numbers.


```python
def problem_2_3_3_solver(n):
    return problem_2_3_1_generic_solver(n, lambda i, j : np.abs(i - j) + 1)

# n = 100
print('n =  100, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_3_solver(100)))
# n = 200
print('n =  200, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_3_solver(200)))
# n = 300
print('n =  300, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_3_solver(300)))
# n = 400
print('n =  400, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_3_solver(400)))
# n = 500
print('n =  500, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_3_solver(500)))
```

    n =  100, forward error = 6.1e-12, error magnification factor = 6.7e+03, condition number = 1e+04
    n =  200, forward error = 3.3e-11, error magnification factor = 1.4e+04, condition number = 4e+04
    n =  300, forward error = 9.9e-11, error magnification factor = 4.4e+04, condition number = 9e+04
    n =  400, forward error = 2e-10, error magnification factor = 8.5e+04, condition number = 1.6e+05
    n =  500, forward error = 2.7e-10, error magnification factor = 1.2e+05, condition number = 2.5e+05


4. Carry out the steps of Computer Problem 3 for the matrix with entries $A_{ij} = \sqrt{(i - j)^2 + n / 10}$.


```python
def problem_2_3_4_solver(n):
    return problem_2_3_1_generic_solver(n, lambda i, j : np.sqrt(np.power(i - j, 2) + n / 10))

# n = 100
print('n =  100, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_4_solver(100)))
# n = 200
print('n =  200, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_4_solver(200)))
# n = 300
print('n =  300, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_4_solver(300)))
# n = 400
print('n =  400, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_4_solver(400)))
# n = 500
print('n =  500, forward error = {:.2g}, error magnification factor = {:.2g}, condition number = {:.2g}'.format(*problem_2_3_4_solver(500)))
```

    n =  100, forward error = 3.7e-09, error magnification factor = 5.1e+06, condition number = 6.2e+07
    n =  200, forward error = 1.3e-06, error magnification factor = 1.4e+09, condition number = 1.3e+10
    n =  300, forward error = 4.6e-05, error magnification factor = 3.5e+10, condition number = 6.2e+11
    n =  400, forward error = 0.00092, error magnification factor = 5.6e+11, condition number = 1.5e+13
    n =  500, forward error = 0.022, error magnification factor = 2.1e+13, condition number = 2.3e+14


5. For what values of n does the solution in Computer Problem 1 have no correct signiﬁcant digits?


```python
print('n = 11, forward error = {:.3g}, error magnification factor = {:.3g}, condition number = {:.3g}'.format(*problem_2_3_1_solver(11)))
```

    n = 11, forward error = 0.0513, error magnification factor = inf, condition number = 4.82e+15


    <ipython-input-28-4e0128b0a498>:30: RuntimeWarning: divide by zero encountered in double_scalars
      error_magnification_factor = relative_forward_error / relative_backward_error


---
# 2.4 The PA=LU Factorization

### Example
Apply Gaussian elimination with partial pivoting to solve the system

\begin{matrix}
  x_1 - x_2 + 3x_3 = -3 & \\ 
-1x_1 - 2x_3 = 1 & \\ 
 2x_1 + 2x_2 + 4x_3 = 0 
\end{matrix}


```python
A = np.array([1, -1, 3, -1, 0, -2, 2, 2, 4]).reshape(3, 3)
b = np.array([-3, 1, 0])
lu, piv = linalg.lu_factor(A)
x = linalg.lu_solve([lu, piv], b)
print(x)
```

    [ 1.  1. -1.]


### Example

Solve the system $2x_1 + 3x_2 = 4$,$3x_1 + 2x_2 = 1$ using the PA = LU factorization with partial pivoting


```python
"""
[[2, 3]
 [3, 2]]
"""
A = np.array([2, 3, 3, 2]).reshape(2, 2)
b = np.array([4, 1])
lu, piv = linalg.lu_factor(A)
x = linalg.lu_solve([lu, piv], b)
print(x)
```

    [-1.  2.]


# 2.5 Iterative Methods

## Jacobi Method


```python
def jacobi_method(A, b, x0, k):
    """
    Use jacobi method to solve equations
    
    Args:
        A  (numpy 2d array): the matrix
        b  (numpy 1d array): the right hand side vector
        x0 (numpy 1d array): initial guess
        k  (real number): iterations
        
    Return:
        The approximate solution
        
    Exceptions:
        ValueError
            The size of matrix's column is not equal to the size of vector's size
    """
    if A.shape[1] is not x0.shape[0] :
        raise ValueError('The size of the columns of matrix A must be equal to the size of the x0')
    
    D = np.diag(A.diagonal())
    inv_D = linalg.inv(D) 
    LU = A - D
    xk = x0
    
    for _ in range(k):
        xk = np.matmul(b - np.matmul(LU, xk), inv_D)
    
    return xk
```

### Example
Apply the Jacobi Method to the system $3u + v = 5$, $u + 2v = 5$


```python
A = np.array([3, 1, 1, 2]).reshape(2, 2)
b = np.array([5, 5])
x = jacobi_method(A, b, np.array([0, 0]), 20)
print('x = %s' %x)
```

    x = [0.99999998 1.99999997]


## Gauss-Seidel Method


```python
def gauss_seidel_method(A, b, x0, k):
    """
    Use gauss seidel method to solve equations
    
    Args:
        A  (numpy 2d array): the matrix
        b  (numpy 1d array): the right hand side vector
        x0 (numpy 1d array): initial guess
        k  (real number): iterations
        
    Return:
        The approximate solution
        
    Exceptions:
        ValueError
            The size of matrix's column is not equal to the size of vector's size
    """
    if A.shape[1] is not x0.shape[0] :
        raise ValueError('The size of the columns of matrix A must be equal to the size of the x0')
    
    D = np.diag(A.diagonal())
    L = np.tril(A) - D
    U = np.triu(A) - D
    inv_LD = linalg.inv(L + D)
    xk = x0
    
    for _ in range(k):
        xk = np.matmul(inv_LD, -np.matmul(U, xk) + b)
    
    return xk
```

### Example
Apply the Gauss-Seidel Method to the system

$$
\begin{bmatrix}
3 & 1 & -1 \\ 
2 & 4 & 1  \\ 
-1 & 2 & 5
\end{bmatrix}
\begin{bmatrix}
u \\ 
v \\ 
w
\end{bmatrix}
=
\begin{bmatrix}
4 \\ 
1 \\ 
1
\end{bmatrix}
$$


```python
A = np.array([3, 1, -1, 2, 4, 1, -1, 2, 5]).reshape(3, 3)
b = np.array([4, 1, 1])
x0 = np.array([0, 0, 0])
gauss_seidel_method(A, b, x0, 24)
```




    array([ 2., -1.,  1.])



## Successive Over-Relaxation


```python
def gauss_seidel_sor_method(A, b, w, x0, k):
    """
    Use gauss seidel method with sor to solve equations
    
    Args:
        A  (numpy 2d array): the matrix
        b  (numpy 1d array): the right hand side vector
        w  (real number): weight
        x0 (numpy 1d array): initial guess
        k  (real number): iterations
        
    Return:
        The approximate solution
        
    Exceptions:
        ValueError
            The size of matrix's column is not equal to the size of vector's size
    """
    if A.shape[1] is not x0.shape[0] :
        raise ValueError('The size of the columns of matrix A must be equal to the size of the x0')
    
    D = np.diag(A.diagonal())
    L = np.tril(A) - D
    U = np.triu(A) - D
    inv_LD = linalg.inv(w * L + D)
    xk = x0
    
    for _ in range(k):
        xk = np.matmul(w * inv_LD, b) + np.matmul(inv_LD, (1 - w) * np.matmul(D, xk) - w * np.matmul(U, xk))
    
    return xk
```

### Example
Apply the Gauss-Seidel Method with sor to the system

$$
\begin{bmatrix}
3 & 1 & -1 \\ 
2 & 4 & 1  \\ 
-1 & 2 & 5
\end{bmatrix}
\begin{bmatrix}
u \\ 
v \\ 
w
\end{bmatrix}
=
\begin{bmatrix}
4 \\ 
1 \\ 
1
\end{bmatrix}
$$


```python
A = np.array([3, 1, -1, 2, 4, 1, -1, 2, 5]).reshape(3, 3)
b = np.array([4, 1, 1])
x0 = np.array([0, 0, 0])
w = 1.25
gauss_seidel_sor_method(A, b, w, x0, 14)
```




    array([ 2., -1.,  1.])



# 2.6 Methods for symmetric positive-definite matrices

##  Cholesky factorization

### Example

Find the Cholesky factorization of 
$\begin{bmatrix}
4 & -2 & 2   \\ 
-2 & 2 & -4  \\ 
2 & -4 & 11
\end{bmatrix}$


```python
A = np.array([4, -2, 2, -2, 2, -4, 2, -4, 11]).reshape(3, 3)
R = linalg.cholesky(A)
print(R)
```

    [[ 2. -1.  1.]
     [ 0.  1. -3.]
     [ 0.  0.  1.]]


## Conjugate Gradient Method


```python
def conjugate_gradient_method(A, b, x0, k):
    """
    Use conjugate gradient to solve linear equations
    
    Args:
        A  : input matrix
        b  : input right hand side vector
        x0 : initial guess
        k  : iteration
        
    Returns:
        approximate solution
    
    
    """
    xk = x0
    dk = rk = b - np.matmul(A, x0)
    for _ in range(k):
        if not np.any(rk) or all( abs(i) <= 1e-16 for i in rk) is True:
            break
        ak = float(np.matmul(rk.T, rk)) / float(np.matmul(dk.T, np.matmul(A, dk)))
        xk = xk + ak * dk
        rk1 = rk - ak * np.matmul(A, dk)
        bk = np.matmul(rk1.T, rk1) / np.matmul(rk.T, rk)
        dk = rk1 + bk * dk
        rk = rk1
    return xk
```

### Example 
Solve 

$$
\begin{bmatrix}
2 & 2 \\ 
2 & 5 \\
\end{bmatrix}
\begin{bmatrix}
u \\ 
v  
\end{bmatrix}
=
\begin{bmatrix}
6 \\ 
3 
\end{bmatrix}
$$

using the Conjugate Gradient Method


```python
A = np.array([2, 2, 2, 5]).reshape(2, 2)
b = np.array([6, 3])
x0 = np.array([0, 0])
conjugate_gradient_method(A, b, x0, 2)
```




    array([ 4., -1.])



### Example 
Solve 

$$
\begin{bmatrix}
1  & -1 & 0 \\
-1 &  2 & 1 \\
0  &  1 & 2 \\
\end{bmatrix}
\begin{bmatrix}
u \\
v \\
w \\
\end{bmatrix}
=
\begin{bmatrix}
0 \\
2 \\
3 \\
\end{bmatrix}
$$


```python
A = np.array([1, -1, 0, -1, 2, 1, 0, 1, 2]).reshape(3, 3)
b = np.array([0, 2, 3])
x0 = np.array([0, 0, 0])
conjugate_gradient_method(A, b, x0, 10)
```




    array([1., 1., 1.])



### Example 
Solve 

$$
\begin{bmatrix}
1  & -1 & 0 \\
-1 &  2 & 1 \\
0  &  1 & 5 \\
\end{bmatrix}
\begin{bmatrix}
u \\
v \\
w \\
\end{bmatrix}
=
\begin{bmatrix}
3  \\
-3 \\
4  \\
\end{bmatrix}
$$


```python
A = np.array([1, -1, 0, -1, 2, 1, 0, 1, 5]).reshape(3, 3)
b = np.array([3, -3, 4])
x0 = np.array([0, 0, 0])
x = slinalg.cg(A, b, x0)[0]
print('x = %s' %x )
```

    x = [ 2. -1.  1.]


## Preconditioning

# 2.7 Nonlinear Systems Of Equations

## Multivariate Newton's Method


```python
def multivariate_newton_method(fA, fDA, x0, k):
    """
    Args:
        fA  (function handle) : coefficient matrix with arguments
        fDA (function handle) : right-hand side vector with arguments
        x0  (numpy 2d array)  : initial guess
        k   (real number)     : iteration
        
    Return:
        Approximate solution xk after k iterations
    """
    xk = x0
    for _ in range(k):
        lu, piv = linalg.lu_factor(fDA(*xk))
        s = linalg.lu_solve([lu, piv], -fA(*xk))
        xk = xk + s
    return xk
```

### Example
Use Newton's method with starting guess $(1,2)$ to find a solution of the system

$$
v - u^3 = 0 \\
u^2 + v^2 - 1 = 0
$$


```python
fA = lambda u,v : np.array([v - pow(u, 3), pow(u, 2) + pow(v, 2) - 1], dtype=np.float64)
fDA = lambda u,v : np.array([-3 * pow(u, 2), 1, 2 * u, 2 * v], dtype=np.float64).reshape(2, 2)
x0 = np.array([1, 2])
multivariate_newton_method(fA, fDA, x0, 10)
```




    array([0.82603136, 0.56362416])



### Example
Use Newton's method to find the solutions of the system

$$
f_1(u,v) = 6u^3 + uv - 3^3 - 4 = 0 \\
f_2(u,v) = u^2 - 18uv^2 + 16v^3 + 1 = 0
$$


```python
fA = lambda u,v : np.array([6 * pow(u, 3) + u * v - 3 * pow(v, 3) - 4,
                           pow(u, 2) - 18 * u * pow(v, 2) + 16 * pow(v, 3) + 1], 
                           dtype=np.float64)

fDA = lambda u,v : np.array([18 * pow(u, 2) + v,  
                            u - 9 * pow(v, 2), 
                            2 * u - 18 * pow(v, 2), 
                            -36 * u * v + 48 * pow(v, 2)], 
                            dtype=np.float64).reshape(2, 2)

x0 = np.array([2, 2], dtype=np.float64)

multivariate_newton_method(fA, fDA, x0, 5)
```




    array([1., 1.])



---
## MIT License

Copyright (c) Jim00000

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

b9f28e29c3ef0d970d3968b9918f4950e88033fc
ipynb
Jupyter Notebook
2_Systems_Of_Equations.ipynb
Jim00000/Numerical-Analysis
b520a9b0e494bea72c3e54e0cae5bbb7c75748c1
MIT
2017-07-08T21:26:07.000Z
2022-02-22T11:09:40.000Z
2_Systems_Of_Equations.ipynb
Jim00000/Numerical-Analysis
b520a9b0e494bea72c3e54e0cae5bbb7c75748c1
MIT
2_Systems_Of_Equations.ipynb
Jim00000/Numerical-Analysis
b520a9b0e494bea72c3e54e0cae5bbb7c75748c1
MIT
2019-01-06T09:55:29.000Z
2021-09-28T15:45:32.000Z
Qwen/Qwen-72B
1. YES
2. YES