diff --git a/functions/convert_matrix_to_triangular_matrix_gauss_naif.m b/functions/convert_matrix_to_triangular_matrix_gauss_naif.m index 9f7ccdd..e48f0f5 100644 --- a/functions/convert_matrix_to_triangular_matrix_gauss_naif.m +++ b/functions/convert_matrix_to_triangular_matrix_gauss_naif.m @@ -15,7 +15,7 @@ function [U,b] = convert_matrix_to_triangular_matrix_gauss_naif (U,b) ok=1; - % convertiamo la matrice in una triangolare superiore + % convertiamo la matrice for i=1:1:n-1 if abs(U(i,i)) <= eps * norma diff --git a/functions/convert_matrix_to_triangular_matrix_gauss_pivoting.m b/functions/convert_matrix_to_triangular_matrix_gauss_pivoting.m index 9bc146f..8c52780 100644 --- a/functions/convert_matrix_to_triangular_matrix_gauss_pivoting.m +++ b/functions/convert_matrix_to_triangular_matrix_gauss_pivoting.m @@ -15,7 +15,7 @@ function [U,b,pivot] = convert_matrix_to_triangular_matrix_gauss_pivoting (U,b) pivot = 1:1:n; pivot = pivot'; - % convertiamo la matrice in una triangolare superiore + % convertiamo la matrice for i=1:1:n-1 x_max = max ( abs(U(i:n,i)) ); diff --git a/functions/linear_system_resolver_lower_triangular_matrix.m b/functions/linear_system_resolver_lower_triangular_matrix.m new file mode 100644 index 0000000..197d555 --- /dev/null +++ b/functions/linear_system_resolver_lower_triangular_matrix.m @@ -0,0 +1,23 @@ +%risoluzione matrice triangolare inferiore con algoritmo +% di sostituzione botton up (all'indietro) + +%example input +% U=[1,0,0;1,1,0;1,5,1]; +% b=[5;-8;-2]; + + +function[x] = linear_system_resolver_lower_triangular_matrix(U,b) + + x=[]; + n= length(b); + x(1) = b(1)/U(1,1); + + for i=2:1:n + somma = 0; + for k=i-1:-1:1 + somma=somma+U(i,k) * x(k); + end + + x(i) = (b(i) - somma)/ U(i,i); + end +end \ No newline at end of file diff --git a/functions/linear_system_resolver_triangular_matrix.m b/functions/linear_system_resolver_upper_triangular_matrix.m similarity index 85% rename from functions/linear_system_resolver_triangular_matrix.m rename to functions/linear_system_resolver_upper_triangular_matrix.m index 65a8613..1be0f7c 100644 --- a/functions/linear_system_resolver_triangular_matrix.m +++ b/functions/linear_system_resolver_upper_triangular_matrix.m @@ -6,7 +6,7 @@ % b=[5;-8;-2]; -function[x] = linear_system_resolver_triangular_matrix(U,b) +function[x] = linear_system_resolver_upper_triangular_matrix(U,b) x=[]; n= length(b); diff --git a/functions/thomas_tridiagonal_matrix_no_pivoting.m b/functions/thomas_tridiagonal_matrix_no_pivoting.m new file mode 100644 index 0000000..afaf7c2 --- /dev/null +++ b/functions/thomas_tridiagonal_matrix_no_pivoting.m @@ -0,0 +1,54 @@ + +% questo algoritmo prende in input una matrice tridiagonale +% e i termini noti : ax=b e ne calcola la soluzione x +% usando il metodo di thomas lineare + +%U=[10,20,30;1,5,3;4,56,34]; matrice di input +%b=[1;2;17]; termini noti + +function [x] = thomas_tridiagonal_matrix_no_pivoting (U,b) + + + n= length(b); % lunghezza del vettore dei termini noti o numero equazioni + norma = norm(U,inf); + + ok=1; + % convertiamo la matrice + for i=1:1:n-1 + + if abs(U(i,i)) <= eps * norma + disp('WARNING: Elemento prossimo allo zero sulla diagonale') + end + + if U(i,i) == 0 + ok=0; + break + + end + j=i+1; + m= U(j,i) / U(i,i); + U(j,i+1) = U(j,i+1) + ( U(i,i+1) * (- m ) ) ; + U(j,i)=m; + + end + + if ok == 0 || U(n,n) == 0 + disp('impossibile convertire') + x=[]; + else + if abs(U(n,n)) <= eps * norma + disp('WARNING: Elemento prossimo allo zero sulla diagonale') + end + + Low=tril(U,-1) + eye (size(U)); + Upp=triu(U); + + %1 ly=b + %2 ux=y + + y = linear_system_resolver_lower_triangular_matrix(Low,b); + x = linear_system_resolver_upper_triangular_matrix(Upp,y); + + end + +end diff --git a/samples/test.m b/samples/sample_matrices.m similarity index 75% rename from samples/test.m rename to samples/sample_matrices.m index 923d553..e31880c 100644 --- a/samples/test.m +++ b/samples/sample_matrices.m @@ -1 +1,5 @@ +%band 2 2 A=[6,-1,-1,0,0,0,0,0,0;-1,6,-1,-2,0,0,0,0,0;-3,-1,6,0,-1,0,0,0,0;0,-1,0,6,-1,-1,0,0,0;0,0,-1,-1,6,-1,-1,0,0;0,0,0,-3,-1,6,0,-1,0;0,0,0,0,-1,0,6,-1,-2;0,0,0,0,0,-1,-1,6,-1;0,0,0,0,0,0,-1,-1,6] + +%tridiagonal +A = [5,1,0,0;4,5,1,0;0,1,3,5;0,0,4,6] diff --git a/samples/spettro.m b/samples/spettro.m new file mode 100644 index 0000000..cf3f9d8 --- /dev/null +++ b/samples/spettro.m @@ -0,0 +1,39 @@ +%Ax=b +%A = X + Y + +%x = inv(X) *b - (inv(X) * Y *x) + +% se C = inv(X) * b e B = -(inv(X) * Y ) +% allora x = Bx + C + +%metodo iterativo converge se e solo se il raggio spettrale di B < 1 +% raggio spettrale = max (abs(eig(B))) + + +%in jacobi X = diag(diag(A)) e Y = A - X +%in gauss-seidel X = tril(A) e Y = A - X + + + +%quindi + +%jacobi B = -(inv(diag(diag(A))) * (A - (diag(diag(A)))) ) +%jacobi C = inv(diag(diag(A))) * b + +%gauss B = -(inv(tril(A)) * (A - (tril(A))) ) +%gauss C = inv(tril(A)) * b + + +%formula definitiva +%jacobi converge se +% max(abs(eig((-(inv(diag(diag(A))) * (A - (diag(diag(A)))) ))))) +% <1 + +%gauss converge se +% max(abs(eig(-(inv(tril(A)) * (A - (tril(A))) )))) +% <1 + + + + +