(* Esercizio 7.3.5 *)

(*a*)
n=3;
a={{3,-1,1},{3,6,2},{3,3,7}};
MatrixForm[a]

a1=a[[1;;1,1;;1]];
a2=a[[1;;2,1;;2]];
a3=a[[1;;3,1;;3]];
MatrixForm[a1]
MatrixForm[a2]
MatrixForm[a3]
Det[a1]
Det[a2]
Det[a3]

b={1,0,4};
x=Table[xx[i],{i,1,n}];
MatrixForm[a];
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

ρ=Max[Abs[λ]]
ω0=2./(1.+Sqrt[1.-ρ^2])

MatrixForm[b]
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=200;tol=1.0 10^(-6);
ω=ω0;
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)+(1.-ω)*xk[k][[i]];
,{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]]

(*c*)
n=4;
a={{10,5,0,0},{5,10,-4,0},{0,-4,8,-1},{0,0,-1,5}};
MatrixForm[a]

a1=a[[1;;1,1;;1]];
a2=a[[1;;2,1;;2]];
a3=a[[1;;3,1;;3]];
MatrixForm[a1]
MatrixForm[a2]
MatrixForm[a3]
Det[a1]
Det[a2]
Det[a3]
Det[a]

b={6,25,-11,-11};
x=Table[xx[i],{i,1,n}];
MatrixForm[a];
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

ρ=Max[Abs[λ]]
ω0=2./(1.+Sqrt[1.-ρ^2])

MatrixForm[b]
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=200;tol=1.0 10^(-6);
ω=ω0-0.01;
ω=0.5;
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)+(1.-ω)*xk[k][[i]];
,{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]]

(* Esercizio 7.3.10 *)

(*a*)
n=3;
a={{1,2,-2},{1,1,1},{2,2,1}};
MatrixForm[a]

a1=a[[1;;1,1;;1]];
a2=a[[1;;2,1;;2]];
a3=a[[1;;3,1;;3]];
MatrixForm[a1]
MatrixForm[a2]
MatrixForm[a3]
Det[a1]
Det[a2]
Det[a3]

b={7,2,5};
x=Table[xx[i],{i,1,n}];
MatrixForm[a];
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

ρJ=Max[Abs[λ]]

(* 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]
λG=Eigenvalues[t]//N

ρG=Max[Abs[λG]]

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=25;tol=1.0 10^(-4);
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]]

(* Esercizio 7.3.14 *)

n=10;
x=Table[xx[i],{i,1,n-1}]

a=Table[0.0,{i,1,n-1},{j,1,n-1}];
b=Table[0.0,{i,1,n-1}];
Do[a[[i,i]]=1.0,{i,1,n-1}];
Do[a[[i,i-1]]=-0.5;a[[i,i+1]]=-0.5,{i,2,n-2}];
a[[1,2]]=a[[n-1,n-2]]=-0.5;
b[[1]]=0.5;

sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[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
ρ=Max[Abs[λ]]
ω0=2./(1.+Sqrt[1.-ρ^2])

Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n-1}];
Do[ff[i]=eq[i][[2]],{i,1,n-1}];
f=Table[ff[i],{i,1,n-1}]

xk[0]=Table[0.0,{i,1,n-1}]
kmax=10;tol=1.0 10^(-4);
k=0;
While[k<kmax,sys1=Table[xx[i]==xk[k][[i]],{i,1,n-1}];
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]]

(* Matrici D, L ed U *)
d=DiagonalMatrix[Diagonal[a]];
Clear[l,u];
l=Table[0.0,{i,1,n-1},{j,1,n-1}];
u=Table[0.0,{i,1,n-1},{j,1,n-1}];
Do[l[[i,j]]=a[[i]][[j]],{i,1,n-1},{j,1,i-1}];
Do[u[[i,j]]=a[[i]][[j]],{i,1,n-1},{j,i+1,n-1}];
t=Dot[Inverse[d-l],u];

MatrixForm[a]
MatrixForm[b]
MatrixForm[d]
MatrixForm[l]
MatrixForm[u]
MatrixForm[t]
λ=Eigenvalues[t]//N
ρ=Max[Abs[λ]]

Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n-1}];
Do[ff[i]=eq[i][[2]],{i,1,n-1}];
f=Table[ff[i],{i,1,n-1}]

xk[0]=Table[0.0,{i,1,n-1}]
kmax=20;tol=1.0 10^(-4);
k=0;
While[k<kmax,
Do[
sysold=Table[xx[j]==xk[k][[j]],{j,i+1,n-1}];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-1}];
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]

Clear[eq,f,ff]
Do[eq[i]=Reduce[sys[[i]],xx[i]],{i,1,n-1}];
Do[ff[i]=eq[i][[2]],{i,1,n-1}];
f=Table[ff[i],{i,1,n-1}]

xk[0]=Table[0.0,{i,1,n-1}]
kmax=20;tol=1.0 10^(-4);
ω=ω0;
k=0;
While[k<kmax,
Do[
sysold=Table[xx[j]==xk[k][[j]],{j,i+1,n-1}];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)+(1.-ω)*xk[k][[i]];
,{i,1,n-1}];
xk[k+1]=Table[xtemp[i],{i,1,n-1}];
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]