For the Ronbrock method, we need to solve a linear system of the form
$$
M_{ij}x_{j}=b_{i} \;,
$$
with M a square matrix (repeated indecies imply summation).

Such systems are soved by (among  other methods) the so-called LU factorization (or decomposition),
where you decompose $M_{ij}=L_{ik}U_{kj}$ with $L_{i, j>i}=0$, $L_{i, j=i}=1$, $U_{i,j<i}=0$.


That is if $M$ is $N \times N$ matrix, L,U are defined as
\begin{align}
    &L=\left( \begin{matrix}
    1 & 0 & 0 & 0 & \dots &0 & 0 \\
    L_{2,1} & 1 & 0 & 0 & \dots &0 & 0\\
    L_{3,1} & L_{3,2} & 1 & 0 & \dots &0 & 0 \\
    \vdots & \vdots & \vdots & \ddots & \dots & \vdots & \vdots \\
    L_{N-1, 1} & L_{N-1, 2} & L_{N-1, 3} & L_{N-1, 4} & \dots & 1 & 0 \\
    L_{N, 1} & L_{N, 2} & L_{N, 3} & L_{N, 4} & \dots & L_{N, N-1} & 1 \\
    \end{matrix}\right) \;, \\
    %
    &U=\left( \begin{matrix}
    U_{1,1} & U_{1,2} & U_{1,3} & U_{1,4} & \dots & U_{1,N-1} & U_{1,N} \\
    0 & U_{2,2} & U_{2,3} & U_{2,4} & \dots & U_{2,N-1} & U_{2,N}\\
    0 & 0 & U_{3,3} & U_{3,4} & \dots & U_{3,N-1} & U_{3,N} \\
    \vdots & \vdots & \vdots & \ddots & \dots & \vdots & \vdots \\
    0 & 0 & 0 &0 & \dots & U_{N-1,N-1} & U_{N-1,N} \\
    0 & 0 & 0& 0 & \dots & 0 &U_{N,N} \\
    \end{matrix}\right)
    %
\end{align}

Then we have in general $M_{i, j} = \sum_{k=1}^{i}L_{i,k}U_{k,j}$. Since 
$L_{i, k \geq i}=0$ and $U_{k>j,j}=0$, the sum runs up to $j$ if $i \geq j$ and
$i$ if $i \leq j$  (for $i=j$ then both are correct). That is




$$
M_{i, j \geq i} = \sum_{k=1}^{i-1}L_{i,k}U_{k,j}+ U_{i,j} \Rightarrow 
U_{i,j }=M_{i,j } - \sum_{k=1}^{i-1}L_{i,k}U_{k,j }\; , \;\;\; j \geq i \\[0.5cm]
M_{i, j \leq i} = \sum_{k=1}^{j-1}L_{i,k}U_{k,j} +L_{i,j}U_{j,j} \Rightarrow 
L_{i,j}=\left( M_{i,j} - \sum_{k=1}^{j-1}L_{i,k}U_{k,j} \right) U_{jj}^{-1} , \;\;\; j \leq i
$$

Since $U$ and $L$ are triangular matrices, we can solve these two systems sequentially
$$
L_{i,k}y_{k}=b_{k} \\
U_{k,j}x_{j}=y_{k},
$$

since 
$$
y_1 = b_{1} \\
L_{2,1}y_{1}+y_{2}=b_{2} \Rightarrow y_{2}=b_{2}-L_{2,1}y_{1} \\
\vdots \\
y_{i}=b_{i} - \sum_{j=1}^{i-1}L_{i,j}y_{j}
$$

and

$$
U_{N,N}x_{N}=y_{N} \Rightarrow x_{N}=y_{N}/U_{N,N}\\
U_{N-1,N}x_{N}+U_{N-1,N-1}x_{N-1}=y_{N-1} \Rightarrow x_{N-1}=\left(y_{N-1} -U_{N-1,N}x_{N} \right)/U_{N-1,N-1} \\
\vdots \\
x_{i}=\dfrac{y_{i} -\displaystyle\sum_{j=i+1}^{N} U_{i,j}x_{j}  }{U_{i,i}}
$$

Since $U_{i,i}$  appears to denominator, if the diagonal terms of $U$ are small (or god forbid they vanish), we would have a problem.


To solve this problem we do $LUP$ decomposition, where $L \; U=P \; M$ with $P$ a permutation matrix so that the diagonal of $U$ has the dominant components in each row.

Then solving $M x =b$ is equavalent to solving $\left( P \; M \right) x =P \; b$ with LU decomposition of 
$P \; M$. That is x solves both systems (no need to permute x).


There is a clever way to make the docomposition faster. This is by initializing 
$L=1_{N \times N}$ and $U=M$. Then we have the follwing algorithm for LU decomposition without
pivoting:


```bash

Input: M, N
#M: matrix
#N: size of M

#initialize U
U=M
#initialize L
L=Unit(N,N)

    for k in [2,...,N] do
        for i in [k,...,N] do
            L[i][k-1]=U[i][k-1]/U[k-1][k-1]
            for j in [k-1,...,N] do
                U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j]
            done
        done
    done
```

I will not write the algorithm including pivoting, as the code in python will not be different.


```python
import numpy as np
```


```python
def Abs(x):
    '''
    Input x  and you'll get the abs.
    '''
    #Not the best way, but it is an easy way to not use numpy
    return (x**2)**0.5

```


```python
def ind_max(row,N):
    '''
    Find the index of the maximum of a list (row) of lentgth N.
    '''
    _in=0
    _max=row[0]
    i=0
    while i<N:#the end of the row should be included (convension in how I use LUP..)
        if row[i]>_max:
            _max=row[i]
            _in=i
        i+=1
            
    return _in
```


```python
def Sum(List,N):
    '''
    Calculates the sum of a List of size N
    '''
    s=0
    for i in range(N):
        s+=List[i]
    return s
```


```python
def index_swap(A,index_1,index_2):
    '''
        index_swap takes an array and interchanges 
         A[index_1] with A[index_2].
         
         Example:
             A=[0,1,2,3]
             index_swap(A,0,2)
             A
             >>[2, 1, 0, 3]
    '''
    
    tmp=A[index_1]
    A[index_1]=A[index_2]
    A[index_2]=tmp
    
    
def apply_permutations_vector(A,P,N):
    '''
    Applies the permutations given by P from LUP
    to a list A of length N, and returns the result.
    Example:
    A=[1,2,5,8,3]
    P=[2,4,0,3,1]

    apply_permutations_vector(A,P,5)
    >>[5, 3, 1, 8, 2]
    '''
    #that is, you make a list like this (P basically gives you the indices of A):)
    
    Ap=[A[ P[i] ] for i in range(N)]

    return Ap
    
def apply_permutations_matrix(M,P,N):
    '''
    Applies the permutations given by P from LUP
    to a N*N array M of length N, and returns the result.
    
    M=[
    [ 0,  2,  2 , 3 , 5],
    [-3, -1,  1 , 5 , 9],
    [ 1, -1,  1 , 4 , 7],
    [ 1, -1,  1 , 0 , 2],
    [ 1, -1,  1 , 0 , 3]
    ]

    P=[2,0,1,4,3]

    apply_permutations_matrix(M,P,5)
    >>[
      [ 1, -1, 1, 4, 7],
      [ 0,  2, 2, 3, 5],
      [-3, -1, 1, 5, 9],
      [ 1, -1, 1, 0, 3],
      [ 1, -1, 1, 0, 2]
      ]
    '''
    Mp=[ [M[ P[i] ][j]for j in range(N)] for i in range(N) ]
    

    return Mp
```


```python
def LUP(M,N,_tiny=1e-20):
    U=[  [ M[i][j] for j in range(N)] for i in range(N) ]
    L=[  [ 0 if i!=j else 1 for j in range(N)] for i in range(N) ]
    #this is the "permutation vector". if it is e.g. [2 1 0 3] it means you make 0<->2
    P=[  i for i in range(N) ]
    
    for k in range(1,N):
        for i in range(k,N):
            #find the index of the maximum in column
            _col=[Abs(U[_r][k-1]) for _r in range(k-1,N)]
            
            #find the index of the maximum of _col
            # notice that the length of _col is N-(k-1)
            len_col=N-(k-1)
            pivot=ind_max( _col ,len_col) + k - 1 #convert the index of _col (it has a length of len_col) to the index of  a row of U   
            
            ##################################################
            #this was in LU_julia (instead of "<_tiny"  it had "== 0").
            #if you remove it, then you get a lot of infinities
            #it has to do with the fact that if U[pivot][k-1] <_tiny , then U[k-1][k-1] will be a zero, 
            #L[i][k-1] explodes. 
            #You are allowed to skip this i, then, because if U[pivot][k-1] <_tiny , then all U[i][k-1] are small!
            #Check that this is true by  uncommenting print(_col)
            if Abs(U[pivot][k-1]) < _tiny  :         
                #print(_col)
                break
            ###################################################
            #if the maximum is not at k-1, swap!
            if pivot != k-1 : 
                # Permute rows k-1 and pivot in U
                
                index_swap(P,k-1,pivot)
                
                tmpU=[U[k-1][_r] for _r in range(k-1,N)]
                
                #print(U)
                for _r in range(k-1,N):
                    U[k-1][_r]=U[pivot][_r]
                #print(U)
                for _r in range(k-1,N):
                    U[pivot][_r]=tmpU[_r-(k-1)]#again we have to convert the index of tmpU
                #print(U)
                #print("=========================")
                tmpL=[L[k-1][_r] for _r in range(k-1)]
                #print(L)
                for _r in range(k-1):
                    L[k-1][_r]=L[pivot][_r]
                #print(L)
                for _r in range(k-1):
                    L[pivot][_r]=tmpL[_r]
                    
                #print(L)
                #print("========================")
                
            L[i][k-1]=U[i][k-1]/U[k-1][k-1]
           
        
            for j in range(k-1,N):
                U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j]

                
    return L,U,P
        
```


```python
def Solve_LU(L,U,P,b,N):
    '''
    This solves P*M*x=P*b (x is also the solution to M*x=b)
    Input:
    L,U,P= LUP decomposition of M. with P*M=L*U
    b=the right hand side of the equation
    N=the number of equations
    '''
    #I know that this uses more memory than what it needs, but I want to keep b as is.
    #(in C++ I will make it a bit beter ;) )
    b=apply_permutations_vector(b,P,N)
    d=[0 for i in range(N) ]
    x=[0 for i in range(N) ]

    d[0]=b[0]
    for i in range(1,N):
        d[i]=b[i]-Sum(  [L[i][j]*d[j] for j in range(i)],i )
    
    
    x[N-1]  = d[N-1]/U[N-1][N-1]  
    for i in range(N-2,-1,-1):
        x[i]=(d[i]-Sum( [U[i][j]*x[j] for j in  range(i+1,N)],N-(i+1) )  )/U[i][i]
        print( Sum( [U[i][j]*x[j] for j in  range(i+1,N)],N-(i+1) )  )
        

    b=apply_permutations_vector(b,P,N)#return b as it was
    return x
```

## tests


```python
#check if Solve_LU works


if True:
    
    NT=5000#NT tests
    N=12#N*N matrices
    testSol=[0 for i in range(NT)]

    for i in range(NT):
        
        #M=np.random.randint(-3,3,size=[N,N])
        b=np.random.rand(N)*13.-6.5
        M=np.random.rand(N,N)*4-2
        L,U,P=LUP(M,N)
        x=Solve_LU(L,U,P,b,N)
        testSol[i]=np.array( np.dot(M,x))-np.array(b)
        
        
        
    print(np.max(testSol))
```

    7.175948724125192e-11



```python
from scipy.linalg import lu_factor,lu_solve,lu

```


```python
#check LUP against numpy.
#in test I will have the maximum difference between my L,U with what np.lu returns, 
#and the difference between my L*U-P*M. So, test should be an array with small numbers!


#even when I get difference with numpy it is not important, because the decomposition is still correct 
#(no nan or inf)!!!!

#change to True to run tests
if True:

    NT=5000#NT tests
    N=12#N*N matrices
    testL=[0 for i in range(NT)]
    testU=[0 for i in range(NT)]
    testM=[0 for i in range(NT)]
    
    
    for i in range(NT):
        
        #M=np.random.randint(-3,3,size=[N,N])
        M=np.random.rand(N,N)*4-2
        L,U,P=LUP(M,N)
        Ps,Ls,Us=lu(M)

        testU[i]=np.max(np.array(U)-Us)
        testL[i]=np.max(np.array(L)-Ls)
        
        testM[i]=np.max(np.dot(L,U)- np.array(apply_permutations_matrix(M,P,N) ))
        if testL[i] > 1e-5:
            #print(np.array(L))
            #print(Ls)
            #print([U[_t][_t] for _t in range(N)])
            print(testM[i])
            pass
            

    print(np.max(testU) , np.max(testL) , np.max(testM))
```

    4.218847493575595e-14 3.7136960173711486e-14 2.220446049250313e-15



```python

```

ae1a857c66eba53405d42d526f78477410296995
ipynb
Jupyter Notebook
Differential_Equations/python/Rosenbrock/jupyter/LU_decomposition.ipynb
dkaramit/ASAP
afade2737b332e7dbf0ea06eb4f31564a478ee40
MIT
Differential_Equations/python/Rosenbrock/jupyter/LU_decomposition.ipynb
dkaramit/ASAP
afade2737b332e7dbf0ea06eb4f31564a478ee40
MIT
Differential_Equations/python/Rosenbrock/jupyter/LU_decomposition.ipynb
dkaramit/ASAP
afade2737b332e7dbf0ea06eb4f31564a478ee40
MIT
2021-12-15T02:03:01.000Z
2021-12-15T02:03:01.000Z
Qwen/Qwen-72B
1. YES
2. YES