% define a matrix >> A = rand(5,5) A = 0.1934 0.6979 0.4966 0.6602 0.7271 0.6822 0.3784 0.8998 0.3420 0.3093 0.3028 0.8600 0.8216 0.2897 0.8385 0.5417 0.8537 0.6449 0.3412 0.5681 0.1509 0.5936 0.8180 0.5341 0.3704 % find its LU decomposition with my program >> [L,U] = LU(A); % check that LU = A >> L*U ans = 0.1934 0.6979 0.4966 0.6602 0.7271 0.6822 0.3784 0.8998 0.3420 0.3093 0.3028 0.8600 0.8216 0.2897 0.8385 0.5417 0.8537 0.6449 0.3412 0.5681 0.1509 0.5936 0.8180 0.5341 0.3704 % find the difference between LU and A >> L*U-A ans = 1.0e-15 * 0 0 0 0 0 0 0 0 0.0555 -0.1110 0 0 0 -0.0555 0 0 -0.1110 0 0.0555 0.1110 0 0 0 0 -0.0555 % find the maximum difference between LU and A >> max(max(L*U-A)) ans = 1.1102e-16 % hunt for matlab's version of the LU decomposition >> help HELP topics: mpugh/matlab - (No table of contents file) matlab/general - General purpose commands. matlab/ops - Operators and special characters. matlab/lang - Programming language constructs. matlab/elmat - Elementary matrices and matrix manipulation. matlab/elfun - Elementary math functions. matlab/specfun - Specialized math functions. matlab/matfun - Matrix functions - numerical linear algebra. matlab/datafun - Data analysis and Fourier transforms. matlab/polyfun - Interpolation and polynomials. matlab/funfun - Function functions and ODE solvers. matlab/sparfun - Sparse matrices. matlab/graph2d - Two dimensional graphs. matlab/graph3d - Three dimensional graphs. matlab/specgraph - Specialized graphs. matlab/graphics - Handle Graphics. matlab/uitools - Graphical user interface tools. matlab/strfun - Character strings. matlab/iofun - File input/output. matlab/timefun - Time and dates. matlab/datatypes - Data types and structures. matlab/demos - Examples and demonstrations. simulink/simulink - Simulink simulink/blocks - Simulink block library. simulink/simdemos - Simulink demonstrations and samples. simulink/dee - Differential Equation Editor toolbox/control - Control System Toolbox. control/obsolete - (No table of contents file) toolbox/local - Preferences. mutools/commands - Mu-Analysis and Synthesis Toolbox. mutools/subs - Mu-Analysis and Synthesis Toolbox. nnet/nnet - Neural Network Toolbox. nnet/nndemos - Neural Network Demonstrations and Applications. toolbox/optim - Optimization Toolbox. toolbox/signal - Signal Processing Toolbox. toolbox/splines - Splines Toolbox. stateflow/stateflow - Stateflow stateflow/sfdemos - (No table of contents file) toolbox/stats - Statistics Toolbox. toolbox/tour - MATLAB Tour For more help on directory/topic, type "help topic". % look in their suggested directory: >> help matlab/matfun Matrix functions - numerical linear algebra. Matrix analysis. norm - Matrix or vector norm. normest - Estimate the matrix 2-norm. rank - Matrix rank. det - Determinant. trace - Sum of diagonal elements. null - Null space. orth - Orthogonalization. rref - Reduced row echelon form. subspace - Angle between two subspaces. Linear equations. \ and / - Linear equation solution; use "help slash". inv - Matrix inverse. cond - Condition number with respect to inversion. condest - 1-norm condition number estimate. chol - Cholesky factorization. cholinc - Incomplete Cholesky factorization. lu - LU factorization. luinc - Incomplete LU factorization. qr - Orthogonal-triangular decomposition. nnls - Non-negative least-squares. pinv - Pseudoinverse. lscov - Least squares with known covariance. Eigenvalues and singular values. eig - Eigenvalues and eigenvectors. svd - Singular value decomposition. eigs - A few eigenvalues. svds - A few singular values. poly - Characteristic polynomial. polyeig - Polynomial eigenvalue problem. condeig - Condition number with respect to eigenvalues. hess - Hessenberg form. qz - QZ factorization for generalized eigenvalues. schur - Schur decomposition. Matrix functions. expm - Matrix exponential. logm - Matrix logarithm. sqrtm - Matrix square root. funm - Evaluate general matrix function. Factorization utilities qrdelete - Delete column from QR factorization. qrinsert - Insert column in QR factorization. rsf2csf - Real block diagonal form to complex diagonal form. cdf2rdf - Complex diagonal form to real block diagonal form. balance - Diagonal scaling to improve eigenvalue accuracy. planerot - Given's plane rotation. % now that I know it's called "lu", I can ask for help by name: >> help lu LU LU factorization. [L,U] = LU(X) stores a upper triangular matrix in U and a "psychologically lower triangular matrix" (i.e. a product of lower triangular and permutation matrices) in L, so that X = L*U. X must be square. [L,U,P] = LU(X) returns lower triangular matrix L, upper triangular matrix U, and permutation matrix P so that P*X = L*U. LU(X), with one output argument, returns the output from LINPACK'S ZGEFA routine. LU(X,THRESH) controls pivoting in sparse matrices, where THRESH is a pivot threshold in [0,1]. Pivoting occurs when the diagonal entry in a column has magnitude less than THRESH times the magnitude of any sub-diagonal entry in that column. THRESH = 0 forces diagonal pivoting. THRESH = 1 is the default. See also LUINC, QR, RREF. % now I use matlab's version of the LU decomposition >> [L,U] = lu(A); % check that LU = A >> L*U ans = 0.1934 0.6979 0.4966 0.6602 0.7271 0.6822 0.3784 0.8998 0.3420 0.3093 0.3028 0.8600 0.8216 0.2897 0.8385 0.5417 0.8537 0.6449 0.3412 0.5681 0.1509 0.5936 0.8180 0.5341 0.3704 % find the difference between LU and A >> L*U-A ans = 1.0e-15 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1110 0 0 % find the maximum difference between LU and A >> max(max(L*U-A)) ans = 1.1102e-16 % Lesson learnt, for this particular matrix, my code and their code % had an equal maximum error, but their code had all zeros at all % other entries while mine didn't. % trying again, with a different matrix: >> A = rand(10) A = Columns 1 through 7 0.4161 0.2417 0.6602 0.9543 0.2919 0.7377 0.6104 0.1864 0.8098 0.3323 0.8814 0.6919 0.0268 0.9886 0.0639 0.9345 0.4768 0.6986 0.4928 0.1022 0.0483 0.0748 0.1288 0.4688 0.3054 0.0835 0.8933 0.9854 0.3100 0.6868 0.7059 0.8293 0.1958 0.7639 0.2047 0.9441 0.2972 0.2399 0.9706 0.9776 0.5827 0.9125 0.9807 0.6472 0.7172 0.3001 0.3662 0.6854 0.6655 0.5551 0.4638 0.8652 0.9981 0.1394 0.9671 0.4623 0.9885 0.9228 0.4110 0.4407 0.0148 0.5503 0.0483 0.6916 0.2417 0.4248 0.0062 0.6406 0.7772 0.4604 Columns 8 through 10 0.8000 0.2033 0.6156 0.2894 0.8193 0.0464 0.6951 0.0584 0.9519 0.2593 0.5385 0.1690 0.7132 0.1902 0.8267 0.7204 0.5995 0.6114 0.7333 0.2923 0.8473 0.6223 0.0913 0.1141 0.9898 0.5068 0.6492 0.1524 0.8841 0.1148 % my version: >> [L,U] = LU(A); % my error: >> max(max(L*U-A)) ans = 3.8303e-15 % matlab's version: >> [L,U] = lu(A); % matlab's error: >> max(max(L*U-A)) ans = 2.2898e-16 >> b = rand(10,1) b = 0.7829 0.0032 0.7970 0.6418 0.1785 0.5294 0.2187 0.5481 0.0582 0.5876 >> [y] = down_solve(L,b) y = 0.7829 -0.3477 1.1214 -0.6360 -1.9253 20.4366 -4.6052 0.3517 4.9958 5.8457 >> [x] = up_solve(U,y) x = -2.0034 -0.1698 0.1666 -1.9730 2.5228 0.9031 0.0437 4.2115 -0.9089 -1.9218 >> max(abs(A*x-b)) ans = 2.8644e-14