Quit[];

Directory[]

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

(* METODO DI JACOBI *)

(* Esercizio 7.3.5 c *)

Clear[a,x]

a=Import["a1.txt","Table"];
b=Import["b1.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

Diagonal[a]

d=Diagonal[a]*IdentityMatrix[n];
r=a-d;
t=-Dot[Inverse[d],r];
c=Dot[Inverse[d],b];
MatrixForm[d]
MatrixForm[r]
MatrixForm[t]
λ=Eigenvalues[t]//N
ρ=Max[Abs[λ]]

kmax=100;tol=1.0 10^(-6);
k=0;
xk[0]=Table[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;
t=-Dot[Inverse[d],r];
c=Dot[Inverse[d],b];
MatrixForm[d]
MatrixForm[r]
MatrixForm[t]
λ=Eigenvalues[t]//N
ρ=Max[Abs[λ]]

kmax=40;tol=1.0 10^(-6);
k=0;xk[0]=Table[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]]