# Exercise


```python
#%pylab inline
%pylab notebook
# pylab Populating the interactive namespace from numpy and matplotlib
# numpy for numerical computation
# matplotlib for ploting
import numpy as np
from numpy import append, array, diagonal, tril, triu
from numpy.linalg import inv
from scipy.linalg import lu
#from scipy.linalg import solve
from pprint import pprint
from numpy import array, zeros, diag, diagflat, dot
from sympy import *
import sympy as sym
init_printing()
```

    Populating the interactive namespace from numpy and matplotlib


    /opt/conda/lib/python3.5/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['Polygon', 'log', 'sqrt', 'source', 'det', 'ones', 'power', 'zeros', 'sinc', 'solve', 'exp', 'interactive', 'product', 'sign', 'tanh', 'sinh', 'plotting', 'test', 'flatten', 'multinomial', 'poly', 'cos', 'Circle', 'cbrt', 'var', 'transpose', 'tan', 'diff', 'roots', 'mod', 'seterr', 'trace', 'plot', 'reshape', 'gamma', 'beta', 'pi', 'eye', 'cosh', 'trunc', 'take', 'diag', 'conjugate', 'floor', 'invert', 'vectorize', 'add', 'prod', 'binomial', 'nan', 'sin']
    `%matplotlib` prevents importing * from pylab and numpy
      "\n`%matplotlib` prevents importing * from pylab and numpy"



```python

```

2.1 Solve Ax =b


```python
A = array([[54, 14, -11, 2],
               [14, 50, -4, 29],
               [-11, -4, 55, 22],
               [2, 29, 22, 95]])
b = array([1, 1, 1, 1])
x0 = [1, 1, 1,1]

```

## (a) L-U decomposition


```python
n = 20

p, l, u = lu(A)
```


```python
p
```




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




```python
l
```




    array([[ 1.        ,  0.        ,  0.        ,  0.        ],
           [ 0.25925926,  1.        ,  0.        ,  0.        ],
           [-0.2037037 , -0.02476038,  1.        ,  0.        ],
           [ 0.03703704,  0.61421725,  0.43831321,  1.        ]])




```python
u
```




    array([[ 54.        ,  14.        , -11.        ,   2.        ],
           [  0.        ,  46.37037037,  -1.14814815,  28.48148148],
           [  0.        ,   0.        ,  52.73083067,  23.11261981],
           [  0.        ,   0.        ,   0.        ,  67.30154198]])




```python
b
```




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




```python
y = np.linalg.solve(l, b)
y    

```




    array([ 1.        ,  0.74074074,  1.22204473, -0.02765113])




```python
x = np.linalg.solve(u, y)
x
```




    array([ 0.01893441,  0.01680508,  0.02335523, -0.00041085])



## (b) Gauss-Jacobi iteration


```python
def gjacobi(A, b, maxit=1000,tol = 1/1000, x0=None):
    """
    Solve a linear equation by the gauss jacobi iteration outlined in the book.
    Follows the eq:
        x = inv(D)(b - Rx)
    Where D is the diagonal matrix of A and R is the remainder s.t D + R = A
    """
    # If we have not provided an initial array for x make a new one
    if x0==None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0
    d = np.diag(np.diag(A))
    
    for it in np.arange(maxit):
        dx =  inv(d).dot(b - A.dot(x))
        x = x + dx
        if np.linalg.norm(dx)<tol:
            break
    #Do not print the result for 1000 times
    #print("x:\n", x) 
    return x

gjacobi(A,b, maxit = 10000, tol = 1/100000)    
```




    array([ 0.01893373,  0.01680349,  0.02335483, -0.00041184])



## (c) Gauss-Seidel iteration
How many Gauss-Jacobi and Gauss-Seidel iterations are required to get answers
that agree with the L-U decomposition solution to four signicant digits?






















```python
def gseidel(A, b, maxit=1000, tol=1/1000,lamb =1, x0=None):
    """
    Solve a linear equation by the gauss seidel iteration outlined in the book
    Follows the eq:
        
    """
    Q = np.tril(A)
    # If we have not provided an initial array for x make a new one
    if x0 == None:
        x = array([1 for _ in range(A.shape[1])])
    else:
        x = x0
    for it in np.arange(maxit):
        dx = np.linalg.inv(Q).dot(b -A.dot(x))
        x = x + lamb*dx
        if np.linalg.norm(dx)<tol:
            break
    #Do not print the result for 1000 times            
    #print("x:\n", x) 
    return x

gseidel(A,b, maxit = 10000, tol = 1/100000)
```




    array([ 0.01893538,  0.0168029 ,  0.02335414, -0.00040996])



### 2.3. 

Use a random number generator create a random 10 by 10 matrix A and a
random 10 by 1 vector b. Then use the Matlab function flop to count the number
of floating point operations needed to solve the linear equation Ax = b 1, 10,
and 50 times for each of the following algorithms:


 - (a) $x = A \ b$
 - (b) $x = U^{-1}(L^{-1}b)$, computing the L-U factors of A only once using the Matlab
function lu.
 - (c) $x = A^{-1}b$, computing $A^{-1}$ only once using function that computes the inverse.

https://docs.python.org/3/library/timeit.html
https://stackoverflow.com/questions/29280470/what-is-timeit-in-python

using IPython’s timeit magic,


```python
%lsmagic
```




    Available line magics:
    %alias  %alias_magic  %autocall  %automagic  %autosave  %bookmark  %cat  %cd  %clear  %colors  %config  %connect_info  %cp  %debug  %dhist  %dirs  %doctest_mode  %ed  %edit  %env  %gui  %hist  %history  %killbgscripts  %ldir  %less  %lf  %lk  %ll  %load  %load_ext  %loadpy  %logoff  %logon  %logstart  %logstate  %logstop  %ls  %lsmagic  %lx  %macro  %magic  %man  %matplotlib  %mkdir  %more  %mv  %notebook  %page  %pastebin  %pdb  %pdef  %pdoc  %pfile  %pinfo  %pinfo2  %popd  %pprint  %precision  %profile  %prun  %psearch  %psource  %pushd  %pwd  %pycat  %pylab  %qtconsole  %quickref  %recall  %rehashx  %reload_ext  %rep  %rerun  %reset  %reset_selective  %rm  %rmdir  %run  %save  %sc  %set_env  %store  %sx  %system  %tb  %time  %timeit  %unalias  %unload_ext  %who  %who_ls  %whos  %xdel  %xmode
    
    Available cell magics:
    %%!  %%HTML  %%SVG  %%bash  %%capture  %%debug  %%file  %%html  %%javascript  %%js  %%latex  %%perl  %%prun  %%pypy  %%python  %%python2  %%python3  %%ruby  %%script  %%sh  %%svg  %%sx  %%system  %%time  %%timeit  %%writefile
    
    Automagic is ON, % prefix IS NOT needed for line magics.




```python
%%timeit
gjacobi(A,b, maxit = 10000, tol = 1/100000)   

```

    1000 loops, best of 3: 575 µs per loop



```python
%%timeit

gseidel(A,b, maxit = 10000, tol = 1/100000)
```

    1000 loops, best of 3: 406 µs per loop



```python

```

f898bb3e732fee55316bf9545493e6a65f23ef94
ipynb
Jupyter Notebook
Exercise02.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
2018-05-09T08:17:44.000Z
2021-12-26T07:02:17.000Z
Exercise02.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
Exercise02.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
2017-12-15T13:39:35.000Z
2021-05-15T15:06:02.000Z
Qwen/Qwen-72B
1. YES
2. YES