MATLAB EXAMPLES Matrix Solution Methods

1. MATLAB EXAMPLESMatrix Solution Methods 58:110 Computer-Aided Engineering Spring 2005 Department of Mechanical and industrial engineering January 2005

2. Some useful functions

3. LU decomposition Here is an example to use LU decomposition for 4X4 matrix A: >> A=[7 3 -1 2;3 8 1 -4; -1 1 4 -1; 2 -4 -1 6] A = 7 3 -1 2 3 8 1 -4 -1 1 4 -1 2 -4 -1 6 >> [l u p]=lu(A) l = 1.0000 0 0 0 0.4286 1.0000 0 0 -0.1429 0.2128 1.0000 0 0.2857 -0.7234 0.0898 1.0000 u = 7.0000 3.0000 -1.0000 2.0000 0 6.7143 1.4286 -4.8571 0 0 3.5532 0.3191 0 0 0 1.8862 p = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 • Note: p is a permutation matrix corresponding to the pivoting strategy used. If the matrix is not diagonally dominant, p will not be an identity matrix.

4. Condition number In MATLAB, Condition number is defined by: To calculate condition number of matrix A by MATLAB built-in function: >> cond(A) ans = 13.7473 Or by the definition, the product of 2-norm of A and inverse A: >> norm(A,2)*norm(inv(A),2) ans = 13.7473

5. Iterative methods – Jocobi method The following M-file shows how to use Jocobi method in MATLAB: (jocobi_example.m) % Iterative Solutions of linear euations:(1) Jocobi Method % Linear system: A x = b % Coefficient matrix A, right-hand side vector b A=[7 3 -1 2; 3 8 1 -4; -1 1 4 -1; 2 -4 -1 6]; b= [-1;0;-3;1]; % Set initial value of x to zero column vector x0=zeros(1,4); % Set Maximum iteration number k_max k_max=1000; % Set the convergence control parameter erp erp=0.0001; % Show the q matrix q=diag(diag(A)) % loop for iterations for k=1:k_max for i=1:4 s=0.0; for j=1:4 if j==i continue else s=s+A(i,j)*x0(j); end end x1(i)=(b(i)-s)/A(i,i); end if norm(x1-x0)<erp break else x0=x1; end end % show the final solution x=x1 % show the total iteration number n_iteration=k

6. Iterative methods- Jocobi method Running M-file in command window: >> jocobi_r_ex q = 7 0 0 0 0 8 0 0 0 0 4 0 0 0 0 6 x = -0.9996 0.9996 -0.9999 0.9996 n_iteration = 60

7. Iterative methods - Gauss-Seidel method In the M-file for Gauss-seidel method, the only difference is the q matrix, which replaced by: q=tril(A), the codes in the loop changed as the following: for i=1:4 s1=0.0; s2=0.0; if (i==1) continue else for j=1:i-1 s1=s1+A(i,j)*x1(j); end end for j=i+1:4 s2=s2+A(i,j)*x0(j); end x1(i)=(b(i)-s1-s2)/A(i,i); end Running M-file in command window >> gauss_example q = 7 0 0 0 3 8 0 0 -1 1 4 0 2 -4 -1 6 x = -0.9996 0.9997 -0.9999 0.9997 n_iteration = 10

8. Iterative methods - SOR method Again, In the M-file for SOR method, the only difference is the q matrix, which replaced by: q1=tril(A)-diag(diag(A)); q2=diag(diag(A))/1.4; q=q1+q2; Note: the relaxation factor set to 1.4 for this case. for i=1:4 s1=0.0; s2=0.0; if (i==1) continue else for j=1:i-1 s1=s1+A(i,j)*x1(j); end end for j=i+1:4 s2=s2+A(i,j)*x0(j); end x1(i)=r*(b(i)-s1-s2)/A(i,i)+(1.0-r)*x0(i); end • Run this M-file: • >> SOR_example q = 5.0000 0 0 0 3.0000 5.7143 0 0 -1.0000 1.0000 2.8571 0 2.0000 -4.0000 -1.0000 4.2857 r = 1.4000 x = -0.9996 0.9997 -0.9999 0.9997 n_iteration = 12

9. Iterative methods – Conjugate gradient method Use cgs, pcg to solve above linear equations, the tolerance is 1.0E-6 (default) >> x=cgs(A,b,0.00001) cgs converged at iteration 4 to a solution with relative residual 1.2e-015 x = -1.0000 1.0000 -1.0000 1.0000 >> x=pcg(A,b,0.00001) pcg converged at iteration 4 to a solution with relative residual 3.1e-016 x = -1.0000 1.0000 -1.0000 1.0000

10. Iterative methods - Summary Note: the exact solution is (-1,1,-1,1)