diff --git a/exercises/test.m b/exercises/test.m new file mode 100644 index 0000000..a93d882 --- /dev/null +++ b/exercises/test.m @@ -0,0 +1,2 @@ + A=[1,3,0;2,5,4;0,5,3] + \ No newline at end of file diff --git a/functions/compact_band_quadratic_matrix.m b/functions/compact_band_quadratic_matrix.m new file mode 100644 index 0000000..2bbc37b --- /dev/null +++ b/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 diff --git a/functions/jacobi_abs_compact.m b/functions/jacobi_abs_compact.m new file mode 100644 index 0000000..4241ea9 --- /dev/null +++ b/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 diff --git a/functions/jacobi_rel_compact.m b/functions/jacobi_rel_compact.m new file mode 100644 index 0000000..3b5eaf9 --- /dev/null +++ b/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 diff --git a/samples/test.m b/samples/test.m index 71a800c..923d553 100644 --- a/samples/test.m +++ b/samples/test.m @@ -1,61 +1 @@ -&% completa risoluzione di un sistema lineare usando la riduzione -% 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 è ') - x -end - +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]