[Content_Types].xml
_rels/.rels
matlab/document.xml
matlab/output.xml
metadata/coreProperties.xml
metadata/mwcoreProperties.xml
metadata/mwcorePropertiesExtension.xml
metadata/mwcorePropertiesReleaseInfo.xml
Iterative Methods: Jacobi Method In this livescript, you will learn how To use the Jacobi method to solve systems of linear equations As with fixed point iteration, our aim is to cast the system [A]\{x\}=\{c\} Equation 1. into the form \{x\}=g(\{x\}) To do this we consider the following decomposition of [A] [A]=[M]-[N] Equation 2. where [M] and [N] are any pair of matrices that satisfy Eq. (2). Therefore, we may write our system as [M]\{x\}=[N]\{x\}+\{c\} Which suggests the iterative formula [M]\{x\}^{(k+1)}=[N]\{x\}^{(k)}+\{c\} This does mean that we have to solve another system of linear equations to determine the \{x\}^{(k+1)} . And if we’re not careful with our choice of [M] , the complexity could be of order \mathcal{O}(kn^3) . However, we have seen that the structure of the linear system can heavily influence the complexity of the algortithm. For example, a tridiagonal matrix solver only requires \mathcal{O}(n) operations to solve (we’ll see this later). So let’s choose [M] to be a matrix that consists of the diagonals of [A] . Then the left-hand-size looks like \pmatrix{a_{11}& & \cr & a_{22} & \cr & & a_{33} \cr & & & \ddots}
\pmatrix{x^{(k+1)}_{1} \cr x^{(k+1)}_{2} \cr x^{(k+1)}_{3} \cr \vdots} And each the next iterate is given by x_{1}^{(k+1)}=\frac{1}{a_{11}}\text{RHS}_{1} x_{2}^{(k+1)}=\frac{1}{a_{22}}\text{RHS}_{2} x_{3}^{(k+1)}=\frac{1}{a_{33}}\text{RHS}_{3} \vdots which is very fast to compute. This gives what is known as the Jacobi method . Part 1: (a) Show that in general, the inverse of a diagonal matrix is computed by replacing each element of the diagonal with its reciprocal. This therefore implies that the \mathcal{O}(n^{3}) numer of operations required by Gaussian elimination may be avoided. We’ll again consider the system of equations 2x-3y+5z=10 4x+7y-2z=-5 2x-4y+25z=31 Equation 3. (b) Show that the matrix [A] is given by [A] = \pmatrix{2 & -3 & 5 \cr 4 & 7 & -2 \cr 2 & -4 & 25} (c) Decompose [A] into [M] and [N] . (d) Show that the iterative formula given by the Jacobi method is \pmatrix{2 & 0 & 0 \cr 0 & 7 & 0 \cr 0 & 0 & 25}
\pmatrix{x^{(k+1)}_{1} \cr x^{(k+1)}_{2} \cr x^{(k+1)}_{3}}
=
–
\pmatrix{0 & -3 & 5 \cr 4 & 0 & -2 \cr 2 & -4 & 0}
\pmatrix{x^{(k)}_{1} \cr x^{(k)}_{2} \cr x^{(k)}_{3}}
+
\pmatrix{10 \cr -5 \cr 31} To extract the diagonal of the matrix [A] , MATLAB has the function \texttt{diag} . If the input is a matrix, then it extracts the diagonal and outputs it as a vector diag([1,2,3; …
4,5,6; …
7,8,9]) However, if the input is a n -vector, \texttt{diag} constructs an n\times n matrix with the input as the diagonal entries. diag([1,2,3]) Therefore the following code may be used to construct the matrices [M] and [N]
. A = [2 -3 5;4 7 -2;2 -4 25];
M = diag(diag(A))
N = M-A We know that the solution is \{x\}=(1,-1,1)^{T} . So let’s consider the initial guess \{x\}_{\text{guess}}=(1,0,0)^{T} (e) Compute \{x\}^{(1)} yourself. Defining the initial guess and \{c\} c = [10;-5;31];
x = [1;0;0]; we may determine \{x\}^{(1)} in MATLAB using x = M^-1*(N*x+c) (f) Run this command a couple more times and convince yourselves that it converges to the true solution \{x\}=(1,-1,1)^{T} . You might have realised that repeatedly pressing run section is quite tiresome and kinda boring. So we should really use a loop to implement the method. As with all iterative methods we’re never sure how many iterations it will take to reach a sufficient accuracy. Therefore it makes most sense to use a while loop until we reach a sufficient accuracy. Define the residual as R(\{x\}^{(k)})=[A]\{x\}^{(k)}-\{c\} which is a measure of how good the solution is. Up until now, we’ve used the \texttt{abs} to measure the size of the residual in the method. However if you apply \texttt{abs} to a vector example_x = [-1 2 4 0 -4];
abs(example_x) we see that \texttt{abs} just takes the absolute value of each array entry. Recall from Linear Algebra the Euclidean distance of a vector \textbf{x} is given by ||\textbf{x}||_{2}=(x_{1}^2+x_{2}^{2}+\cdots)^{\frac{1}{2}} This is what is called a 2-norm, which may easily be generalised into a p-norm by ||\textbf{x}||_{p}=(x_{1}^{p}+x_{2}^{p}+\cdots)^{\frac{1}{p}} which may be called in MATLAB using \texttt{norm(x,p)} , where \texttt{x} is the vector and \texttt{p} is the norm. NOTE: The \infty -norm is given by ||\textbf{x}||_{\infty}=\max_{i}{|x_{i}|} Therefore, using the 2-norm, the Jacobi method is given by A = [2 -3 5;4 7 -2;2 -4 25];
c = [10;-5;31];
M = diag(diag(A))
N = M-A
x = [1;0;0];
i = 1;
tol = 10^-6;
fprintf(‘ x1 x2 x3 \n’)
fprintf(‘%10.7f %10.7f %10.7f\n’,x(1),x(2),x(3))
while norm(A*x-c,2)>tol
x = M^-1*(N*x+c);
fprintf(‘%10.7f %10.7f %10.7f\n’,x(1),x(2),x(3))
i = i+1;
end
fprintf(‘Iteration was completed in %d iterations’,i) Just as a note, you shouldn’t use this code in an assignment as I have just gotten MATLAB to invert the matrix [M] , and you should never do that whenever there is a better method.
manual code ready 0.4 true false diag([1,2,3; …
4,5,6; …
7,8,9]) 0 55 57 false true diag([1,2,3]) 1 59 59 true false A = [2 -3 5;4 7 -2;2 -4 25]; 2 62 62 false false M = diag(diag(A)) 3 63 63 false true N = M-A 4 64 64 true false c = [10;-5;31]; 5 68 68 false true x = [1;0;0]; 6 68 68 true true x = M^-1*(N*x+c) 7 74 74 true false example_x = [-1 2 4 0 -4]; 8 83 83 false false abs(example_x) 9 84 84 false false A = [2 -3 5;4 7 -2;2 -4 25]; 10 94 94 false false c = [10;-5;31]; 11 95 95 false false M = diag(diag(A)) 12 97 97 false false N = M-A 13 98 98 false false x = [1;0;0]; 14 100 100 false false i = 1; 15 101 101 false false tol = 10^-6; 16 102 102 false false fprintf(‘ x1 x2 x3 \n’) 17 104 104 false false fprintf(‘%10.7f %10.7f %10.7f\n’,x(1),x(2),x(3)) 18 105 105 false false while norm(A*x-c,2)>tol
x = M^-1*(N*x+c);
fprintf(‘%10.7f %10.7f %10.7f\n’,x(1),x(2),x(3))
i = i+1;
end 19 107 111 false true fprintf(‘Iteration was completed in %d iterations’,i) 20 113 113
2020-03-19T09:19:00Z 2020-03-20T06:17:17Z
application/vnd.mathworks.matlab.code MATLAB Code R2019b
9.7.0.1183724 de1efc99-f0d7-49b7-9a57-9a5ce29afbe6
9.7.0.1296695
R2019b
Update 4
Jan 20 2020
3546228467