Numerical Methods  

LUP decompositionThe stability of LU decomposition is improved if pivoting is used to maximize the absolute values of the diagonal elements of the upper triangular matrix $U$. The function LUP_decomp(A) performs LUdecomposition with partial pivoting. It is the same as the function LU_matrix(A) except a pivoting section has been added. It returns an object consisting of the LU matrix, the permutation matrix, and the number of row exchanges made during partial pivoting. It is clear that with the strategy of partial pivoting we avoid the problem caused by the presence of a zero element on the diagonal of $U$. The only exception occurs if during the calculation all the $U_{ij},\quad i=j,\ldots,n$ are simultaneously zero. In that case the matrix of coefficients is singular and the system of equations cannot be solved with this method. Matrix inversionThe inverse of matrix A is defined by the equation, \begin{equation} A \cdot A^{1} = I, \end{equation}where $I$ is the identity matrix. The column vectors of $A^{1}$ can be found by solving $A\cdot\vec{x}=\vec{b}$ where $\vec{b}$ is the corresponding column vector of the identity matrix. We use LUPdecomposition to find the permutation matrix $P$ and the upper and lower triangular matricies $L$ and $U$ such that $P\cdot A= L\cdot U$. Then we can call LU_solve($LU$,$P\cdot I$) to find $A^{1}$. function inverse(A) { // returns the inverse of a square matrix LUP = LUP_decomp(A); Pb = dot(LUP.P,identity(A.length)); Ainv = LU_solve(LUP.LU,Pb); return Ainv; } Calculation of the determinant of a matrixThe determinant of the product of two matrices is the product of their determinants. After LUPdecomposition we have, $\text{det}(P)\text{det}(A)=\text{det}(L)\text{det}(U)$. The determinant of the permutation matrix $P$ is $1^n$ where $n$ is the number of row exchanges made during partial pivoting. The determinant of an upper or lower triangular matrix is the product of its diagonal elements: $\text{det}(L)=1$, $\text{det}(U)=\prod\limits_{i=1}^NU_{ii}$. The following code will calculate the determinant. function det(A) { //returns the determinant of square matrix A LUP = LUP_decomp(A); D = 1; for (i=0; i < A.length; i++) {//product the diagonal elements of the LU matrix D = D*LUP.LU[i][i]; } D = D*LUP.sgn; return D; }
