jacobi and gauss-seidel

Hi everyone, can you say where I am wrong please?

convergence with Jacobi and Gauss-Seidel

A=[p+1 1/2 0 0; 1 p+1 1/3 0; 0 1/2 p+1 1/4; 0 0 1/3 p+1]

p=1

A=[p+1 1/2 0 0; 1 p+1 1/3 0; 0 1/2 p+1 1/4; 0 0 1/3 p+1]

p=[1]

alpha=ones(4,1);

b=A*alpha;

omega=1

[n,m]=size(A);

D=diag(diag(A));

BJ=eye(n)-omega*inv(D)*A;

rhoJ=max(abs(eig(BJ)))

D=diag(diag(A));

OE=omega*tril(A,-1);

BGS=eye(n)-omega*inv(D+OE)*A;

rhoG=max(abs(eig(BGS)))

x0=[1 0 0 0]'; nmax=30;toll=1e-9;

[x,iter,residuo,rho]=Gauss_Seidel_ril(A,b,x0,omega,nmax,toll);

for i=1:iter+1

err(i)=norm(x(i,:)-alpha',inf);

end

it=[0:iter]';

tab=[it err' res ];

fprintf('iter errore residuo\n');

fprintf('%3d %10.2f %10.2f %10.2f %10.2f %10.2e\n')

[John D'Errico]

You should understand this won't work:

A=[p+1 1/2 0 0; 1 p+1 1/3 0; 0 1/2 p+1 1/4; 0 0 1/3 p+1]
p=[1]

p does not exist when you created A. Therefore unless you have defined p before, then the code you wrote will throw an error.

[John D'Errico]

Need I also point out the fallacy in using a Gauss-Seidel code that calls the function inv? Lol. Be serious. Do you think there is not a fundamental flaw in what you are doing, if you write this?

BGS=eye(n)-omega*inv(D+OE)*A; 

If you are going to use inv, you can completely avoid it all.

x = inv(A)*b;

What are the odds that won't be acceptable either? Better yet would be to write it as

x = A\b;

but we won't get into that.

[Chiara]

Thank you.

But BGS=eye(n)-omega*inv(D+OE)*A is a iteration matrix instead x = A\b is a vector, in Bgs inv is used to construct the Gauss-Seidel matrix

[John D'Errico]

No. You do not understand. Using inv to perform Gauss-Seidel is a major problem. Since you could have used inv to solve the entire problem, while never needing to use Gauss-Seidel at all, then using inv in the middle of the algorithm is simply incorrect. It invalidates the entire algorithm you are trying to code.