Quit[];

(* METODO DI JACOBI *)

(* Esempio 1. *)

Clear[a,x]
n=4;
a={{10,5,0,0},{5,10,-4,0},{0,-4,8,-1},{0,0,-1,5}};
b={6,25,-11,-11};
x=Table[xx[i],{i,1,n}];
Det[a];

(* Matrici D ed R *)
d=DiagonalMatrix[Diagonal[a]];
r=a-d;
t=-Dot[Inverse[d],r];
MatrixForm[a]
MatrixForm[b]
MatrixForm[d];
MatrixForm[r];
MatrixForm[t];
Eigenvalues[t]//N

sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[1]]//N
Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n}];
Do[ff[i]=eq[i][[2]],{i,1,n}];
f=Table[ff[i],{i,1,n}];
xk[0]=Table[0.0,{i,1,n}];
kmax=100;tol=1.0 10^(-4);
k=0;
While[k<kmax,sys1=Table[xx[i]==xk[k][[i]],{i,1,n}];
sol1=Solve[sys1,x][[1]];xk[k+1]=f/.sol1//N;
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]

(* Esempio 2. *)
Clear[a,x]
n=3;
a={{1,5,0},{5,-1,-4},{0,-4,0.5}};
b={6,25,-11};
x=Table[xx[i],{i,1,n}];

(* Matrici D ed R *)
d=DiagonalMatrix[Diagonal[a]];
r=a-d;
t=-Dot[Inverse[d],r];
MatrixForm[a]
MatrixForm[b]
MatrixForm[d];
MatrixForm[r];
MatrixForm[t];
Eigenvalues[t]//N

sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[1]]//N
Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n}];
Do[ff[i]=eq[i][[2]],{i,1,n}];
f=Table[ff[i],{i,1,n}];
xk[0]={-2.,1.6,-9.1};
xk[0]=Table[0.0,{i,1,n}];
kmax=10;tol=1.0 10^(-6);
k=0;
While[k<kmax,sys1=Table[xx[i]==xk[k][[i]],{i,1,n}];
sol1=Solve[sys1,x][[1]];xk[k+1]=f/.sol1//N;
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]

(* METODO DI GAUSS-SEIDEL *)

(* Esempio 1. *)
Clear[a,x]
n=4;
a={{10,5,0,1},{5,10,-4,0},{0,-4,8,-1},{-1,0,-1,5}};
b={6,25,-11,-11};
x=Table[xx[i],{i,1,n}];
Det[a]

(* Matrici D, L ed U *)
d=DiagonalMatrix[Diagonal[a]]

Clear[l,u];
l=Table[0.0,{i,1,n},{j,1,n}];
u=Table[0.0,{i,1,n},{j,1,n}];
Do[l[[i,j]]=a[[i]][[j]],{i,1,n},{j,1,i-1}];
Do[u[[i,j]]=a[[i]][[j]],{i,1,n},{j,i+1,n}];
t=Dot[Inverse[d-l],u];

MatrixForm[a]
MatrixForm[b]
MatrixForm[d];
MatrixForm[l];
MatrixForm[u];
MatrixForm[t];
Eigenvalues[t]//N

sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[1]]//N
Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n}];
Do[ff[i]=eq[i][[2]],{i,1,n}];
f=Table[ff[i],{i,1,n}];
xk[0]=Table[0.0,{i,1,n}];
kmax=100;tol=1.0 10^(-6);
k=0;
While[k<kmax,
Do[
sysold=Table[xx[j]==xk[k][[j]],{j,i+1,n}];sysnew=Table[xx[j]==xtemp[j],{j,1,i-1}];
sys1=Join[sysnew,sysold];
sol1=Solve[sys1,x][[1]];
(*Print[i,sysold,sysnew,sys1,sol1];*)
xtemp[i]=f[[i]]/.sol1//N;
,{i,1,n}];
xk[k+1]=Table[xtemp[i],{i,1,n}];
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]

(* Esempio 2. *)
Clear[a,x]
n=3;
a={{1,5,0},{5,-1,-4},{0,-4,0.5}};
b={6,25,-11};
x=Table[xx[i],{i,1,n}];

(* Matrici D, L ed U *)
d=DiagonalMatrix[Diagonal[a]]

Clear[l,u];
l=Table[0.0,{i,1,n},{j,1,n}];
u=Table[0.0,{i,1,n},{j,1,n}];
Do[l[[i,j]]=a[[i]][[j]],{i,1,n},{j,1,i-1}];
Do[u[[i,j]]=a[[i]][[j]],{i,1,n},{j,i+1,n}];
t=Dot[Inverse[d-l],u];

MatrixForm[a]
MatrixForm[b]
MatrixForm[d];
MatrixForm[l];
MatrixForm[u];
MatrixForm[t];
Eigenvalues[t]//N

sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[1]]//N
Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n}];
Do[ff[i]=eq[i][[2]],{i,1,n}];
f=Table[ff[i],{i,1,n}];
xk[0]=Table[0.0,{i,1,n}];
kmax=10;tol=1.0 10^(-6);
k=0;
While[k<kmax,
Do[
sysold=Table[xx[j]==xk[k][[j]],{j,i+1,n}];sysnew=Table[xx[j]==xtemp[j],{j,1,i-1}];
sys1=Join[sysnew,sysold];
sol1=Solve[sys1,x][[1]];
(*Print[i,sysold,sysnew,sys1,sol1];*)
xtemp[i]=f[[i]]/.sol1//N;
,{i,1,n}];
xk[k+1]=Table[xtemp[i],{i,1,n}];
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]