thomas tridiagonal and spectral radius

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

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)) );
 

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

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);

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

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]

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
+
+
+
+
+