question: Prove that if \(A\) has this property (unit row diagonal dominant)
\[ a_{ii}=1>{\displaystyle \sum \limits _{\substack {j=1\\j\neq i}}^{n}} \left \vert a_{ij}\right \vert \ \ \ \ \ \ \ \ \left (1\leq i\leq n\right ) \]
then Richardson iteration is successful.
Solution:
Since the iterative formula is
\begin {align*} x_{k+1} & =x_{k}+Q^{-1}\left (b-Ax_{k}\right ) \\ & =x_{k}+Q^{-1}b-Q^{-1}Ax_{k}\\ & =(I-Q^{-1}A)x_{k}+Q^{-1}b \end {align*}
This converges, by theorem (1) on page 210 when \(\Vert I-Q^{-1}A\Vert <1\)
In Richardson method, \(Q=I\), hence Richardson converges if \(\Vert I-A\Vert <1\)
But \(\Vert I-A\Vert =\Vert \begin {bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ 0 & \cdots & \ddots & \vdots \\ 0 & \cdots & \cdots & 1 \end {bmatrix}\begin {bmatrix} 1 & a_{12} & \cdots & a_{1n}\\ a_{21} & 1 & \cdots & a_{2n}\\ \vdots & \cdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & 1 \end {bmatrix} \Vert =\Vert \begin {bmatrix} 0 & a_{12} & \cdots & a_{1n}\\ a_{21} & 0 & \cdots & a_{2n}\\ \vdots & \cdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & 0 \end {bmatrix} \Vert \)
But since row unit diagonal dominant, then the sum of each row elements remaining (after \(a_{ii}\) was annihilated) is a sum which is less than 1. Hence each row about will sum to some value which is less than 1. Hence the infinity norm of the above matrix, which is the maximum row sum, is less than 1. Hence
\[ \Vert I-A\Vert <1 \]
Hence Richardson will converge. Each iteration will move closer to the solution \(x^{\ast }\)
Problem: Prove that the eigenvalues of a Hermitian matrix are real.
Answer: \(A\) is Hermitian if \(\overline {\left (A^{T}\right ) }=A\,\ \), where the bar above indicates taking the complex conjugate. Hence the matrix is transposed and then each element will be complex conjugated.
Now, an eigenvalue of a matrix is defined such as\[ Ax=\lambda x \]
pre mutliply both sides by \(\overline {\left (x^{T}\right ) }\)
\begin {align*} \overline {\left (x^{T}\right ) }Ax & =\overline {\left (x^{T}\right ) }\lambda x\\ \overline {\left (x^{T}\right ) }Ax & =\lambda \overline {\left ( x^{T}\right ) }x \end {align*}
But since \(A\) is Herminitian, then \(\overline {\left (A^{T}\right ) }=A\), hence the above becomes
\begin {align*} \overline {\left (x^{T}\right ) }\overline {\left (A^{T}\right ) }x & =\lambda \overline {\left (x^{T}\right ) }x\\ \overline {\left (x^{T}A^{T}\right ) }x & =\lambda \overline {\left ( x^{T}\right ) }x\\ \overline {\left (Ax\right ) ^{T}}x & =\lambda \overline {\left ( x^{T}\right ) }x \end {align*}
But \(Ax=\lambda x\), hence the above becomes
\begin {align*} \overline {\left (\lambda x\right ) ^{T}}x & =\lambda \overline {\left ( x^{T}\right ) }x\\ \overline {\left (x^{T}\lambda \right ) }x & =\lambda \overline {\left ( x^{T}\right ) }x\\ \overline {\left (x^{T}\right ) }\overline {\lambda }x & =\lambda \overline {\left (x^{T}\right ) }x\\ \overline {\lambda }\overline {\left (x^{T}\right ) }x & =\lambda \overline {\left (x^{T}\right ) }x\\ \overline {\lambda } & =\lambda \end {align*}
Hence since complex conjugate of eigenvalue is the same as the eigenvalue, therefor \(\lambda \) is real.
Problem: Prove that if \(A\) is nonsingular, then \(A\overline {A^{T}}\) is positive definite.
Answer: \(A\) is nonsingular, meaning its left and right inverses exist and are the same. To show that a matrix is positive definite, need to show that \(x^{T}Ax>0\) for all \(x\neq 0.\)
Let \(\overline {A^{T}}=B\), then let \(n=x^{T}A\overline {A^{T}}x\) \(,\) we need to show that \(n>0\)\[ n=x^{T}\left (AB\right ) x \]
Since \(A^{-1}\) exist, then multiply both sides by \(A^{-1}\) and \(B^{-1}\)
\begin {align*} A^{-1}B^{-1}n & =A^{-1}B^{-1}x^{T}\left (AB\right ) x\\ & =x^{T}x \end {align*}
But \(x^{T}x=\Vert x\Vert ^{2}>0\) unless \(x=0\), hence above becomes
\begin {align*} A^{-1}B^{-1}n & >0\\ A^{-1}\left (\overline {A^{T}}\right ) ^{-1}n & >0\\ \left (\overline {A^{T}}A\right ) ^{-1}n & >0 \end {align*}
Multiply both sides by \(\overline {A^{T}}A\) leads to \(n>0\), hence \(A\overline {A^{T}}\) is positive definite.
Problem: Prove that if \(A\) is positive definite, then its eigenvalues are positive\(.\)
Answer: \(A\) is positive definite implies \(x^{T}Ax>0\) for all \(x\neq 0\). i.e. \(\left \langle x,Ax\right \rangle >0\).
But \(Ax=\lambda x\,\), hence \[ \left \langle x,\lambda x\right \rangle >0 \] or \[ \lambda \left \langle x,x\right \rangle >0 \] But \(\left \langle x,x\right \rangle =\Vert x\Vert ^{2}>0\) unless \(x=0\), therefor \[ \lambda >0 \]
question: Prove that if \(A\) is symmetric, then the gradient of the function \(q\relax (x) =\left \langle x,Ax\right \rangle -2\left \langle x,b\right \rangle \) at \(x\) is \(2\left (Ax-b\right ) \). Recall that the gradient of a function \(g:\mathbb {R} ^{n}\rightarrow \mathbb {R} \) is the vector whose components are \(\frac {\partial g}{\partial x_{i}}\) for \(i=1,2,\cdots n\)
answer:
\(\left \langle a,b\right \rangle \) is \(a^{T}b\), hence using this definition, we can expend the RHS above and see that it will give the result required.
\begin {align*} \left \langle x,Ax\right \rangle & =x^{T}\left (Ax\right ) \\ & =\left [ x_{1}\ x_{2}\ \cdots \ x_{n}\right ] \begin {bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn}\end {bmatrix}\begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} \\ & =\left [ x_{1}\ x_{2}\ \cdots \ x_{n}\right ] \begin {bmatrix} a_{11}x_{1}+a_{12}x_{2}+\cdots +a_{1n}x_{n}\\ a_{21}x_{1}+a_{22}x_{2}+\cdots +a_{2n}x_{n}\\ \vdots \\ a_{n1}x_{1}+a_{n2}x_{2}+\cdots +a_{nn}x_{n}\end {bmatrix} \\ & =x_{1}\left (a_{11}x_{1}+a_{12}x_{2}+\cdots +a_{1n}x_{n}\right ) +x_{2}\left (a_{21}x_{1}+a_{22}x_{2}+\cdots +a_{2n}x_{n}\right ) \\ & +x_{3}\left (a_{31}x_{1}+a_{32}x_{2}+\cdots +a_{3n}x_{n}\right ) +\cdots +x_{n}\left (a_{n1}x_{1}+a_{n2}x_{2}+\cdots +a_{nn}x_{n}\right ) \\ & =\left (a_{11}x_{1}^{2}+a_{12}x_{1}x_{2}+\cdots +a_{1n}x_{1}x_{n}\right ) +\left (a_{21}x_{2}x_{1}+a_{22}x_{2}^{2}+\cdots +a_{2n}x_{2}x_{n}\right ) \\ & +\left (a_{31}x_{1}x_{3}+a_{32}x_{2}x_{3}+a_{33}x_{3}^{2}\cdots +a_{3n}x_{n}x_{3}\right ) +\cdots +\left (a_{n1}x_{n}x_{1}+a_{n2}x_{n}x_{2}+\cdots +a_{nn}x_{n}^{2}\right ) \end {align*}
And
\begin {align*} \left \langle x,b\right \rangle & =x^{T}b\\ & =\left [ x_{1}\ x_{2}\ \cdots \ x_{n}\right ] \begin {bmatrix} b_{1}\\ b_{2}\\ \vdots \\ b_{n}\end {bmatrix} \\ & =x_{1}b_{1}+x_{2}b_{2}+\ \cdots +x_{n}b_{n} \end {align*}
Hence Putting the above together, we obtain
\begin {align*} q\relax (x) & =\left \langle x,Ax\right \rangle -2\left \langle x,b\right \rangle \\ & =\left (a_{11}x_{1}^{2}+a_{12}x_{1}x_{2}+\cdots +a_{1n}x_{1}x_{n}\right ) +\left (a_{21}x_{2}x_{1}+a_{22}x_{2}^{2}+\cdots +a_{2n}x_{2}x_{n}\right ) \\ & +\left (a_{31}x_{1}x_{3}+a_{32}x_{2}x_{3}+a_{33}x_{3}^{2}+\cdots +a_{3n}x_{n}x_{3}\right ) +\cdots +\left (a_{n1}x_{n}x_{1}+a_{n2}x_{n}x_{2}+\cdots +a_{nn}x_{n}^{2}\right ) \\ & -2\left (x_{1}b_{1}+x_{2}b_{2}+\ \cdots +x_{n}b_{n}\right ) \end {align*}
Now taking the derivative of the above w.r.t \(x_{1},x_{2},\cdots ,x_{n}\) to generate the gradient vector, we obtain
\begin {align*} \frac {\partial q\relax (x) }{\partial x_{1}} & =\left (2a_{11}x_{1}+a_{12}x_{2}+\cdots +a_{1n}x_{n}\right ) +\left (a_{21}x_{2}\right ) +\left (a_{31}x_{3}\right ) +\cdots +\left (a_{n1}x_{n}\right ) -2\left ( b_{1}\right ) \\ \frac {\partial q\relax (x) }{\partial x_{2}} & =\left (a_{12}x_{1}\right ) +\left (a_{21}x_{1}+2a_{22}x_{2}+\cdots +a_{2n}x_{n}\right ) +\left (a_{32}x_{3}\right ) +\cdots +\left (a_{n2}x_{n}\right ) -2\left ( b_{2}\right ) \\ & \vdots \\ \frac {\partial q\relax (x) }{\partial x_{n}} & =\left (a_{1n}x_{1}\right ) +\left (a_{2n}x_{2}\right ) +\left (a_{3n}x_{3}\right ) +\cdots +\left (a_{n1}x_{1}+a_{n2}x_{2}+\cdots +2a_{nn}x_{n}\right ) -2\left ( b_{n}\right ) \end {align*}
Hence
\begin {align*} \frac {\partial q\relax (x) }{\partial x_{1}} & =\left [ 2a_{11},\ a_{12,}\ \cdots ,a_{1n}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} +\left [ 0\ ,a_{21},a_{31},\ \cdots ,a_{n1}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{1}\right ) \\ \frac {\partial q\relax (x) }{\partial x_{2}} & =\left [ a_{21}\ ,2a_{22}\ ,\cdots ,a_{2n}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} +\left [ a_{12}\ ,0,\ \cdots ,a_{n2}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{2}\right ) \\ & \vdots \\ \frac {\partial q\relax (x) }{\partial x_{n}} & =\left [ a_{n1}\ ,a_{n2},\ \cdots ,2a_{nn}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} +\left [ a_{n1},\ a_{n2,}\ \cdots ,0\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{n}\right ) \end {align*}
Combine, we get
\begin {align*} \frac {\partial q\relax (x) }{\partial x_{1}} & =\left [ 2a_{11},\ 2a_{12,}\ \cdots ,2a_{1n}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{1}\right ) \\ \frac {\partial q\relax (x) }{\partial x_{2}} & =\left [ 2a_{21}\ ,2a_{22}\ ,\cdots ,2a_{2n}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{2}\right ) \\ & \vdots \\ \frac {\partial q\relax (x) }{\partial x_{n}} & =\left [ 2a_{n1}\ ,2a_{n2},\ \cdots ,2a_{nn}\right ] \begin {bmatrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end {bmatrix} -2\left (b_{n}\right ) \end {align*}
Hence, in Vector/Matrix notation we obtain
\begin {align*} \frac {\partial q\relax (x) }{\partial \vec {x}} & =2A\vec {x}-2\vec {b}\\ & =2\left (A\vec {x}-\vec {b}\right ) \end {align*}
Which is what we are asked to show.
(it would also have been possible to solve the above by expressing everything in the summation notations, i.e. writing \(\left (Ax\right ) _{i}=\sum _{j}^{n}A\left (i,j\right ) x_{j}\), and apply differentiations directly on these summation expression as is without expanding them as I did above, it would have been probably shorter solution)
Question: Prove that the minimum value of \(q\relax (x) \) is \(-\left \langle b,A^{-1}b\right \rangle \)
Solution: \[ q\relax (x) =\left \langle x,Ax\right \rangle -2\left \langle x,b\right \rangle \]
From problem (1) we found that \(\frac {\partial q\relax (x) }{\partial \vec {x}}=2\left (A\vec {x}-\vec {b}\right ) \), hence setting this to zero
\begin {align*} A\vec {x}-\vec {b} & =0\\ \vec {x} & =A^{-1}b \end {align*}
To check if this is a min,max, or saddle, we differentiate \(\frac {\partial q\relax (x) }{\partial \vec {x}}\) once more, and plug in this solution and check
\begin {align*} \frac {\partial }{\partial \vec {x}}\left (\frac {\partial q\relax (x) }{\partial \vec {x}}\right ) & =2\frac {\partial }{\partial \vec {x}}\left ( A\vec {x}-\vec {b}\right ) \\ & =2A \end {align*}
(\(A\) needs to be positive definite here, problem did not say this?) Hence this is a minimum.
Hence minimum value \(q\relax (x) \) is \[ q_{\min }\relax (x) =\left \langle A^{-1}b,Ax\right \rangle -2\left \langle A^{-1}b,b\right \rangle \]
But \(Ax=b\,\), hence
\begin {align*} q_{\min }\relax (x) & =\left \langle A^{-1}b,b\right \rangle -2\left \langle A^{-1}b,b\right \rangle \\ & =-\left \langle A^{-1}b,b\right \rangle \\ & =-\left \langle b,A^{-1}b\right \rangle \end {align*}
question: Prove that if \(\hat {t}=\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }\),and if \(y=x+\hat {t}v\) then \(\left \langle v,b-Ay\right \rangle =0\)
answer: \[ y=x+\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }v \]
Then \begin {align*} \left \langle v,b-Ay\right \rangle & =\left \langle v,b-A\left (x+\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }v\right ) \right \rangle \\ & =\left \langle v,b-Ax-A\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }v\right \rangle \\ & =\left \langle v,b-Ax-Av\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }\right \rangle \\ & =v^{T}\left (b-Ax-Av\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }\right ) \\ & =v^{T}b-v^{T}Ax-v^{T}Av\frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }\\ & =v^{T}b-v^{T}Ax-\left \langle v,Av\right \rangle \frac {\left \langle v,b-Ax\right \rangle }{\left \langle v,Av\right \rangle }\\ & =v^{T}b-v^{T}Ax-\left \langle v,b-Ax\right \rangle \\ & =v^{T}\left (b-Ax\right ) -\left \langle v,b-Ax\right \rangle \\ & =\left \langle v,b-Ax\right \rangle -\left \langle v,b-Ax\right \rangle \\ & =0 \end {align*}
Jacobi, Gauss-Seidel, and SOR iterative solvers were implemented (in Matlab) and compared for rate of convergence. The implementation was tested on the following system
\[ A=\begin {bmatrix} -4 & 2 & 1 & 0 & 0\\ 1 & -4 & 1 & 1 & 0\\ 2 & 1 & -4 & 1 & 2\\ 0 & 1 & 1 & -4 & 1\\ 0 & 0 & 1 & 2 & 04 \end {bmatrix}\begin {bmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\\ x_{5}\end {bmatrix} =\begin {bmatrix} -4\\ 11\\ -16\\ 11\\ -4 \end {bmatrix} \]
and the result was compared to Matlab \(A\backslash b^{\prime }\) result, which is
A\b' ans = 1.00000000000000 -2.00000000000000 4.00000000000000 -2.00000000000000 1.00000000000000 >>
These iterative solvers solve the linear system \(Ax=b\) by iteratively approaching a solution. In all of these methods, the general iteration formula is
\[ x_{k+1}=x_{k}+Q^{-1}\left (I-Ax_{k}\right ) \]
Each method uses a different \(Q\) matrix. For Jacobi, \(Q=diag(A)\), and for Gauss-Seidel, \(Q=L\relax (A) \), which is the lower triangular part of \(A\). For SOR, \(Q=\frac {1}{\omega }\left (diag\relax (A) +\omega L^{0}\relax (A) \right ) \) where \(L^{0}\relax (A) \) is the strictly lower triangular part of \(A\), and \(\omega \) is an input parameter generally \(0<\omega <2\).
For the Jacobi and Gauss-Seidel methods, we are guaranteed to converge to a solution if \(A\) is diagonally dominant. For SOR the condition of convergence is \(\rho \relax (G) <1\), where \(G=\left (I-Q^{-1}A\right ) \), and \(\rho \relax (G) \) is the spectral radius of \(G\).
The following table shows the result of running Jacobian and Gauss-Seidel on the same Matrix. The columns are the relative error between each successive iteration. Defined as \(\frac {\left \vert x_{k+1}-x_{k}\right \vert }{\left \vert x_{k+1}\right \vert }\). The first column is the result from running the Jacobian method, and the second is the result from the Gauss-Seidel method.
| Iteration | Jacobi | Gauss-Seidel |
| \(1\) | \(0.94197873843414\) | \(0.93435567469174\) |
| \(2\) | \(0.22011159808466\) | \(0.13203098317070\) |
| \(3\) | \(0.04613055928361\) | \(0.03056375635161\) |
| \(4\) | \(0.02582530967034\) | \(0.02878251571226\) |
| \(5\) | \(0.02166103846735\) | \(0.02294426703511\) |
| \(6\) | \(0.01508738881859\) | \(0.01935275191913\) |
| \(7\) | \(0.01447950687767\) | \(0.01618762392092\) |
| \(8\) | \(0.01246529285252\) | \(0.01353611057259\) |
| \(9\) | \(0.01166876946105\) | \(0.01130586303057\) |
| \(10\) | \(0.01055706870353\) | \(0.00943545181644\) |
| \(11\) | \(0.00971262454554\) | \(0.00786920716886\) |
| \(12\) | \(0.00886458421973\) | \(0.00655940732143\) |
| \(13\) | \(0.00811604581680\) | \(0.00546521547469\) |
| \(14\) | \(0.00741643318343\) | \(0.00455191661070\) |
| \(15\) | \(0.00677994422374\) | \(0.00379012917894\) |
| \(16\) | \(0.00619437598238\) | \(0.00315507291046\) |
| \(17\) | \(0.00565882948442\) | \(0.00262590555197\) |
| \(18\) | \(0.00516810821828\) | \(0.00218513512307\) |
| \(19\) | \(0.00471915312779\) | \(0.00181810678260\) |
| \(20\) | \(0.00430837834153\) | \(0.00151255968550\) |
We see from the above table that Gauss-Seidel method has faster convergence than Jacobi method. The following is a plot of the above table.
For the SOR method, it depends on the value of \(\omega \) used. First we look at SOR on its own, comparing its performance as \(\omega \) changes. Next, we pick 2 values for \(\omega \) and compare SOR using these values with the Jacobian and the Gauss-Seidel method.
This table below shows the values of the relative error for difference \(\omega \). Only the first 20 iterations are shown. The following \(\omega \) values are used: \(.25,.5,.75,1,1.25,1.5,1.75\). This table also shows the number of iterations needed to achieve the same error tolerance specified by the user. Smaller number of iterations needed to reach this error limit indicates that \(\omega \) selected was better than otherwise.
From the above table, and for the specific data used in this exercise, we observe that \(\omega =1.25\) achieved the best convergence
This is a plot of the above table. (using zoom to make the important part more visible).
Now we show the difference between Jacobi, Gauss-Seidel and SOR (using \(\omega =1.25\)) as this \(\omega \) gave the best convergence using SOR. The following is a plot showing that SOR with \(\omega =1.25\) achieved best convergence.
The following is an output from one run generated by a Matlab script written to test these implementations. 2 test cases used. A,b input as given in the HW assignment sheet given in the class, and a second A,b input shown in the textbook (SPD matrix, but not diagonally dominant) shown on page 245. This is the output.
Conclusion SOR with \(\omega =1.25\) achieved the best convergence. However, one needs to find determine which \(\omega \) can do the best job, and this also will depend on the input. For different input different \(\omega \) might give better result. Hence to use SOR one must first do a number of trials to determine the best value to use. Between the Jacobian and the Gauss-Seidel methods, the Gauss-Seidel method showed better convergence.
In these operations long count, since these algorithms are iterative, hence the operation count will be based on one iteration. One would then need to multiply this count by the number of iterations need to converge as per the specification that the user supplies. Gauss-Seidel will require less iterations to converge to the same solution as compared to Jacobi method. With the SOR method, it will depend on the \(\omega \) selected.
Jacobi Long operations count for one iteration of the Jacobi method.
while keepLookingForSolution k=k+1; xnew=xold+Qinv*(b-A*xold);- - - - - (1) currentError=norm(xnew-xold);- - - - - (2) relError(k)=currentError/norm(xnew);- - - - - (3) if norm(b-A*xnew)<=resLimit || currentError<=errorLimit || k>maxIter keepLookingForSolution=FALSE; else xold=xnew; end end
where \[ Qinv=eye(nRow)/diag(diag(A)); \]
For line(1): For multiplying an \(n\times n\) matrix by \(n\times 1\) vector, \(n^{2}\) ops are required. Hence for \(A\ast xold\) we require \(n^{2}\) ops. The result of \((b-A\ast xold)\) is an \(n\) long vector. Next we have \(Qinv\) multiplied by this \(n\) long vector. This will require only \(n\) ops since \(Qinv\) is all zeros except at the diagonal. Hence to calculate xnew we require \(n+n^{2}\) ops.
For line(2,3): It involves 2 norm operations and one division. Hence we need to find the count for the norm operation. Assuming norm 2 is being used which is defined as \(\sqrt {\sum _{i=1}^{n}x_{i}^{2}}\) hence this requires \(n\) multiplications and one long operation added for taking square root (ok to do?). Hence the total ops for finding relative error is \(2n+2+1=2n+3\)
The next operation is the check for the result limit:
\[ norm(b-A\ast xnew)<=resLimit \]
This involves \(n^{2}\) operations for the matrix by vector multiplication, and \(n+1\) operations for the norm. Hence \(n^{2}+n+1\)
Adding all the above we obtain: \(n+n^{2}+2n+3+n^{2}+n+1=2n^{2}+4n+4\) hence this is \(O\left (2n^{2}\right ) \)
The above is the cost per one iteration. Hence if \(M\) is close to \(n\), then it is worst than G.E. But if the number of iteration is much less than \(N\), then this method will be faster than non-iterative methods based on Gaussian elimination which required \(O\left (n^{3}\right ) \).
Gauss-Seidel Long operations count for one iteration of the Gauss-Seidel. For the statement
while keepLookingForSolution k=k+1; xnew=xold; for i=1:nRow xnew(i)=xnew(i)+ (b(i)-A(i,:)*xnew)/A(i,i);- - - - -(1) end currentError=norm(xnew-xold); - - - - -(2) relError(k)=currentError/norm(xnew); - -- - - -(3) if norm(b-A*xnew)<=resLimit || currentError<=errorLimit || k>maxIter keepLookingForSolution=FALSE; else xold=xnew; end end
For line (1): For each i, there is \(n\) for the operation A(i,:)*xnew, and one operation for the division. Hence the total cost is \(n^{2}\) since there are \(n\) rows.
For line(2+3): The cost is the same as with the Jacobi method, which was found to be \(2n+3\)
The next operation is the check for the result limit:
\[ norm(b-A\ast xnew)<=resLimit \]
This involves \(n^{2}\) operations for the matrix by vector multiplication, and \(n+1\) operations for the norm.
Hence the total operation is \(n^{2}+2n+3+n+1=\)\(n^{2}+3n+4\) hence this is \(O\left (n^{2}\right ) \) per one iteration.
SOR This is the same as Gauss-Seidel, except there is an extra multiplication by \(\omega \) per one row per one iteration. Hence we need to add \(n\) to the count found in Gauss-Seidel. Hence the count for SOR is \(n^{2}+4n+4\) or \(O\left (n^{2}\right ) \)
The following is the output from testing the Steepest descent iterative algorithm.
while keepLookingForSolution k=k+1; v=b-A*xold; - - - - - - (1) t=dot(v,v)/dot(v,A*v); - - - - -(2) xnew=xold+t*v; - - - - -(3) currentError=norm(xnew-xold); - - - - -(4) relError(k)=currentError/norm(xnew); - - - -(5) if norm(b-A*xnew)<=resLimit || currentError<=errorLimit || k>maxIter keepLookingForSolution=FALSE; else xold=xnew; end end
For line (1): \(n^{2}\)
For line (2): For numerator: dot operation is \(n.\) For denominator: for A*v need \(n^{2}\), and then need \(n\) added to the dot operation, hence need \(n+n^{2}+n\) operations. Add 1 to the division, hence need \(n^{2}+2n+1\) for line (2).
For line(3): \(n\) multiplications.
For line(4): \(n+1\) for the norm.
For line(5): \(n+2\) operations.
And finally for norm(b-A*xnew) in the final check, requires \(n^{2}+n+1\)
Hence total is \(n^{2}+n^{2}+2n+1+n+n+1+n+2+n^{2}+n+1=3n^{2}+6n+5\)