Quit[];

SetDirectory["/Users/Lucio/Didattica/Analisi Numerica/Mathematica files"]

(* METODO SOR *)

(* Esercizio 7.3.5 c *)

Clear[a,x]
a=Import["a1.txt","Table"];
b=Import["b1.txt","Table"][[1]];
(*a=Import["a5c.txt","Table"];
b=Import["b5c.txt","Table"][[1]];*)
n=Length[b];
x=Table[xx[i],{i,1,n}];
MatrixForm[a]
Det[a]
MatrixForm[b]
sys=Thread[Dot[a,x]==b]
xsol=x/.Solve[sys,x][[1]]//N

(* E' a definita positiva ? *)

a1=a[[1;;1,1;;1]]
a2=a[[1;;2,1;;2]]
MatrixForm[a1]
MatrixForm[a2]

Do[s[i]=a[[1;;i,1;;i]];Print[Det[s[i]]],{i,1,n}]

MatrixForm[a]
MatrixForm[s[2]]

d=Diagonal[a]*IdentityMatrix[n];
r=a-d;
l=u=Table[0.,{i,1,n},{j,1,n}];
Do[l[[i,j]]=-a[[i,j]],{i,2,n},{j,1,i-1}]
Do[u[[i,j]]=-a[[i,j]],{i,1,n-1},{j,i+1,n}]
MatrixForm[d]
MatrixForm[l]
MatrixForm[u]
tj=-Dot[Inverse[d],r];
λj=Eigenvalues[tj]//N;
ρj=Max[Abs[λj]]

ω0=2.0/(1.0+Sqrt[1.0-ρj^2])

ω=ω0;
t=Dot[Inverse[d-ω*l],ω*u+(1-ω)*d]
c=Dot[ω*Inverse[d-ω*l],b];
MatrixForm[d];
MatrixForm[l];
MatrixForm[u];
MatrixForm[t]
MatrixForm[c]
λ=Eigenvalues[t]//N;
ρ=Max[Abs[λ]]

kmax=20;tol=1.0 10^(-6);
k=0;xk[0]=Table[0.0,{i,1,n}];
While[k<kmax,
xk[k+1]=Dot[t,xk[k]]+c;
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]

(* Esempio 2. *)

Clear[a,x]
a=Import["a.txt","Table"];
b=Import["b.txt","Table"][[1]];
n=Length[b];
x=Table[xx[i],{i,1,n}];
MatrixForm[a]
Det[a]
MatrixForm[b]
sys=Thread[Dot[a,x]==b]
xsol=x/.Solve[sys,x][[1]]//N

d=Diagonal[a]*IdentityMatrix[n];
r=a-d;
l=u=Table[0,{i,1,n},{j,1,n}];
Do[l[[i,j]]=-a[[i,j]],{i,2,n},{j,1,i-1}]
Do[u[[i,j]]=-a[[i,j]],{i,1,n-1},{j,i+1,n}]
MatrixForm[d]
MatrixForm[l]
MatrixForm[u]
tj=-Dot[Inverse[d],r];
λj=Eigenvalues[tj]//N;
ρj=Max[Abs[λj]]
ω0=2.0/(1.0+Sqrt[1.0-ρj^2])

ω=1.1;
t=Dot[Inverse[d-ω*l],ω*u+(1-ω)*d]
c=Dot[ω*Inverse[d-ω*l],b];
MatrixForm[d];
MatrixForm[l];
MatrixForm[u];
MatrixForm[t]
MatrixForm[c]
λ=Eigenvalues[t]//N;
ρ=Max[Abs[λ]]

kmax=20;tol=1.0 10^(-6);
k=0;xk[0]=Table[0.0,{i,1,n}]
While[k<kmax,
xk[k+1]=Dot[t,xk[k]]+c;
norm=Max[Abs[xk[k+1]-xk[k]]];Print[k, "  ", norm];If[norm≤tol,Break[]];
k++]
Print[xk[k]]