compattazione di matrice a banda e jacobi

exercises/test.m

@@ -0,0 +1,2 @@
+ A=[1,3,0;2,5,4;0,5,3]
+ 

functions/compact_band_quadratic_matrix.m

@@ -0,0 +1,43 @@
+%conversione matrice quadrata a banda in versione compatta
+
+function [A] = compact_band_quadratic_matrix (U,left,right)
+
+dim = size(U);
+dimy= dim(2);
+A = zeros(left + right + 1,dimy);
+
+%upper
+x=1;Ax=right;
+for y=2:(right + 1)
+    
+    Uy=y;
+    for Ux=x:(dimy - (y - 1))
+        
+        A(Ax,Uy) = U(Ux,Uy);
+      
+        Uy=Uy+1;
+    end
+    
+    
+     Ax=Ax-1; 
+end
+
+%lower
+y=1;Ax=right + 2 ;
+for x = 2:(left +1)
+    
+    Ux=x;
+    for Uy=y: (dimy - (x- 1))
+       A(Ax,Uy)=U(Ux,Uy); 
+      
+       Ux=Ux+1;
+    end
+    
+   Ax=Ax+1; 
+end
+
+%diagonal
+for x=1:dimy
+    A(right+1,x) = U(x,x);
+    
+end

functions/jacobi_abs_compact.m

@@ -0,0 +1,68 @@
+
+% jacobi algorithm
+
+%input
+%U= matrice di input compattata
+%b=[1;2;17]; termini noti
+%x0 vettore di partenza
+%toll tolleranza assoluta , esempio 1e-4
+%nmax numero massimo di passaggi, esempio 500
+
+
+%Ooutput
+%x0 soluzione computata
+%err errore assoluto computato
+%niter numero di interazioni eseguite
+%ier 1 se il numero di iterazioni ha raggiunto nmax, 0 altrimenti
+
+function [x0,err,niter,ier] = jacobi_abs_compact (U,b,x0,toll,nmax)
+
+    n= length(b); % lunghezza del vettore dei termini noti o numero equazioni
+    niter = 0;
+    ier = 0;
+    err = inf;% err infinito per superare subito toll
+    x1 = x0;% solo per inizializzare
+    
+    %detect right and left
+    x=1;
+    y=1;
+    left= 0;
+    right= 0;
+    while U(x,y)==0
+        right = right + 1;
+        x=x+1;
+    end
+    dim= size(U);
+    dimx=dim(1);
+    left = dimx  - right - 1;
+    
+    %begin iterations
+    
+    while ( niter < nmax ) && ( err >= toll)
+        for i=1:n
+            partial_sum = 0;
+            partial_sum2 = 0;
+
+            kend = min (left, i - 1);
+            for k=1:kend
+                partial_sum = partial_sum + ( U(right +1 +k,i-k) * x0(i-k));
+            end
+            
+            kend = min (right, n - i);
+            for k=1: kend
+                partial_sum2 = partial_sum2 + ( U(right+1-k,i+k) * x0(i+k));
+            end
+            x1(i) = (b(i) - partial_sum - partial_sum2)/U(right +1,i) ;
+        end
+        niter = niter + 1;
+        err = norm(x1 - x0,inf);
+        x0=x1;
+        mex =  [' Iterazione ', num2str(niter),' : ',mat2str(x0), ' Errore assoluto : ',  num2str(err)];
+        disp (mex)
+    end
+    if niter == nmax
+        disp('Warning: Massimo numero di step raggiunti')
+        ier = 1;
+    end
+ 
+end

functions/jacobi_rel_compact.m

@@ -0,0 +1,68 @@
+
+% jacobi algorithm
+
+%input
+%U= matrice di input compattata
+%b=[1;2;17]; termini noti
+%x0 vettore di partenza
+%toll tolleranza relativa , esempio 1e-4
+%nmax numero massimo di passaggi, esempio 500
+
+
+%Ooutput
+%x0 soluzione computata
+%err errore relativo computato
+%niter numero di interazioni eseguite
+%ier 1 se il numero di iterazioni ha raggiunto nmax, 0 altrimenti
+
+function [x0,err,niter,ier] = jacobi_rel_compact (U,b,x0,toll,nmax)
+
+    n= length(b); % lunghezza del vettore dei termini noti o numero equazioni
+    niter = 0;
+    ier = 0;
+    err = inf;% err infinito per superare subito toll
+    x1 = x0;% solo per inizializzare
+    
+    %detect right and left
+    x=1;
+    y=1;
+    left= 0;
+    right= 0;
+    while U(x,y)==0
+        right = right + 1;
+        x=x+1;
+    end
+    dim= size(U);
+    dimx=dim(1);
+    left = dimx  - right - 1;
+    
+    %begin iterations
+    
+    while ( niter < nmax ) && ( err >= toll)
+        for i=1:n
+            partial_sum = 0;
+            partial_sum2 = 0;
+
+            kend = min (left, i - 1);
+            for k=1:kend
+                partial_sum = partial_sum + ( U(right +1 +k,i-k) * x0(i-k));
+            end
+            
+            kend = min (right, n - i);
+            for k=1: kend
+                partial_sum2 = partial_sum2 + ( U(right+1-k,i+k) * x0(i+k));
+            end
+            x1(i) = (b(i) - partial_sum - partial_sum2)/U(right +1,i) ;
+        end
+        niter = niter + 1;
+        err = norm(x1 - x0,inf) / norm (x1,inf);
+        x0=x1;
+        mex =  [' Iterazione ', num2str(niter),' : ',mat2str(x0), ' Errore relativo : ',  num2str(err)];
+        disp (mex)
+    end
+    if niter == nmax
+        disp('Warning: Massimo numero di step raggiunti')
+        ier = 1;
+    end
+ 
+end

samples/test.m

@@ -1,61 +1 @@
-&% completa risoluzione di un sistema lineare usando la riduzione
+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]
-% in matrice triangolare superiore con gauss pivoting parziale
-% e algoritmo di bottom up con sostituzione
-
-U=[0,20,30;0,5,3;0,56,34]; % matrice di input
-b=[1;2;17]; %termini noti
-
-
-
-
-n= length(b); 
-U= [U,b]; 
-U
-
-
-for i=1:1:n-1
-    x_max = max ( abs(U(i:n,i)) );
-    if x_max == 0
-        disp('errore matrice di input, det 0')
-        break
-    else
-        [x,y]= ind2sub(size(U),  find (abs(U(i:n,i)) == x_max) );
-        x= x + i -1;
-        y = i;
-        
-       if x~= i
-            U([i x],:) = U([x i],:);
-            
-            
-       end
-        
-       for j=i+1:1:n
-
-            U(j,:) =  U(j,:) + ( U(i,:) * (- U(j,i) / U(i,i) ) ) ; 
-
-       end
-    end
-end
-
-if x_max ~= 0
-
-    b=U(:,n+1);
-    U=U(1:n,1:n);
-
-
-    x=[];
-    n= length(b);
-    x(n) = b(n)/U(n,n);
-
-    for i=n-1:-1:1
-        somma = 0;
-        for k=i+1:n
-            somma=somma+U(i,k) * x(k);
-        end
-
-        x(i) = (b(i) - somma)/ U(i,i);
-    end
-    disp ('La soluzione <EFBFBD> ')
-    x
-end
-