# Matlab project - Jacobi method for tridiagonal matrices…

I have to do a project in Matlab to my University and I don't quite understand what I should do. I was given script that solves systems of equations with Jacobi's method with given tolerance and number of iterations and I was told to use it.

The title of the project is:

"Solving system of linear equations $Cz=c$ where $C=A+iB$, $A, B\in \mathbb{R}^{n\times n}$ are tridiagonal matrices, $c=a+ib$, $z=x+iy$, $\ x,y,a,b\in\mathbb{R}^n$ with Jacobi method applied to a matrix $$M=\begin{pmatrix} A&-B \\ B&A \end{pmatrix}\in\mathbb{R}^{2n\times 2n}$$

There should be only three diagonals of $A$ and $B$ kept in the computer's memory (in the form of vectors). Compare results with Matlab build-in function."

I don't quite understand what do I need M matrix for. I was given a hint that:

$Cz=c \Leftrightarrow (A+iB)(x+iy)=a+ib \Leftrightarrow (Ax-By)+i(Bx+Ay)=a+ib \Leftrightarrow \begin{pmatrix} A & -B\\ B & A \end{pmatrix}\begin{pmatrix} x\\ y \end{pmatrix} = \begin{pmatrix} a\\ b \end{pmatrix} \Leftrightarrow$

If we introduce following denotation: $\begin{pmatrix} x\\ y \end{pmatrix}=z_1$ and $\begin{pmatrix} a\\ y \end{pmatrix}=c_1$

we can end with one more equivalence:

$\Leftrightarrow Mz_1=c_1$

But now I don't know where did the imaginary unit go and I if need to include it in my code or not.

Here is my code for Jacobi method (x0 is a vector of initial guesses, tol is tolerance, max_iter is maximal interation, we solve $Ax=b$):

function [ x,k ] = Jacobi2( A,b,x0,tol,max_iter )

poprawka=tol+1;
k=0; D=diag(diag(A)); % czy detA!=0?
while(norm(poprawka)>tol) && k<=max_iter
poprawka=D\(b-A*x0);
x=x0+poprawka;
k=k+1;
x0=x;
end
end


Now what I'm doing is that I take those three diagonals, compute matrix $M$ and with given $x,y$ I am trying to solve $Mz=\begin{pmatrix} x\\ y \end{pmatrix}$ ($z$ is intitial guess vector) with above function but I don't think I am getting the right results.

Second thing is that I want to write a function generating random, convergent (in the sense of Jacobi Method, that means $\max(eig(D^{-1})<1$ where $D^{-1}$ is inverse diagonal of matrix $M$ and "eig" are eigenvalues. I know how to generate random matrix $M$ but how to generate convergent matrix M?

Here is my code for generating random $M, A, B$ matrices and vectors $x, y$ (I call them $z_1, z_2$ in my function), every number is in the interval $[0,40]$:

function [ M, z1, z2, A, B ] = projekt_random_tridiagonal( n )
$n - dimension of A and B x=40; A=diag(diag(x*rand(n,n)))+diag((x*rand(n-1,1)),1)+diag((x*rand(n-1,1)),-1); B=diag(diag(x*rand(n,n)))+diag((x*rand(n-1,1)),1)+diag((x*rand(n-1,1)),-1); M=[A -B; B A]; z1=x*rand(n,1); z2=x*rand(n,1); end  I will appreciate any kind of help? The main problem is understanding the question (tried asking my professor and just got the hint above). • To solve this:$(Ax−By)+i(Bx+Ay)=a+ib$you need that the real and imaginary parts are equal, so you have to solve$Ax-By=a$and$Bx+Ay=b$simultaneously. Also, you should be able to replace D\(b-A*x0) with (b-A*x0)./diag(A). – David Dec 1 '14 at 3:45 • @David, could you tell me how to solve for example$Ax-By=a$using the code I was given? I can't think of any way to implement that into my programme. All I can think of is transforming that equation into$x=A^{-1}a+A^{-1}By$but I don't know what to do next since there are two unknowns. Also, what about the$M$matrix and applying it into my code? I feel really stupid right now, but we have never operated on complex matrices before, neither in lectures, nor in labs so it confuses me a lot. – Mateusz Dec 1 '14 at 9:23 • Can you see how the two equations are equivalent to that matrix multiplication$Mz_1=c_1$? You don't really have to think about them as complex matrices once you have split the original equation into 2. – David Dec 1 '14 at 9:31 • Ah, I see it. So I have to solve this system now:$ \left\{\begin{matrix} Ax-By=a\\ Bx+Ay=b \end{matrix}\right \Leftrightarrow Mz_1=c_1$, right? So in my programme I should ask the user to input three diagonals of$A,B$, then I should ask for$x,y,a,b$and I should compute then matrix$M$and vectors$\begin{pmatrix} x\\y \end{pmatrix}=z_1$and$\begin{pmatrix} a\\b \end{pmatrix}=c_1\$ and then use my Jacobi2 function? – Mateusz Dec 1 '14 at 9:44
• Yeah that's right. – David Dec 1 '14 at 10:50