[MAtlab] help su 2 esercizi

[gandelf]

Salve a tutti, ho difficoltà ad implementare 2 esercizi:

1) Scrivi un algoritmo in Matlab che risolva il sistema lineare Ax=b in cui A è una matrice a banda (p,q) determinando prima le matrici U ed L con l’algoritmo su indicato e poi risolvendo Ly=b e Ux=y.
Sperimenta l’algoritmo con n=60, q=p=0, n/4, n/2, n-1. Calcola il tempo di esecuzione per ciascuno di questi parametri mediante la funzione di MatLab cputime.
Suggerimento: Il problema a caso può essere determinato, generando prima una matrice A, casuale di dimensione n i cui elementi sono tra 0 e 1, poi utilizzando successivamente le funzioni predefinite tril(A,k) e triu(A,k) (l’help di triu() e tril() indicano quale valore dare a k). Si sceglie, ugualmente a caso, il vettore colonna xx di dimensione n e quindi si calcola b=A*xx. n, A e b sono i dati del problema; xx è la soluzione esatta del problema.

L’algoritmo di Gauss per la fattorizzazione A=L*U di matrici a banda inferiore p e superiore q è:

for k=1,2,…,n-1
for i=k+1,…,min(k+p,n)
A(i,k) = A(i,k) / A(k,k)
end
for i=k+1,…,min(k+p,n)
for j=k+1,…,min(k+q,n)
A(i,j) = A(i,j) - A(i,k)*A(k,j)
end
end
end


2)Risolvi un sistema lineare Ax=b utilizzando sia il metodo iterativo di Jacobi sia il metodo di Gauss-Seidel.

Suggerimento. I due metodi sono iterativi; indicando con k l’indice di iterazione si ha che la convergenza è ottenuta per k che tende all’infinito. Pertanto nella formulazione dell’algoritmo è necessario inserire un criterio di stop. Questo può essere: fissata un’accuratezza acc termina le iterazioni quando la norma della differenza di due n-ple successive x(k+1), x(k) è minore di acc. Inoltre poiché il criterio di stop potrebbe non essere mai soddisfatto ( caso della non convergenza), entrando in un ciclo infinito, conviene aggiungere un altro criterio di stop: l’esecuzione dell’algoritmo è interrotta quando viene raggiunta un’iterazione massima prefissata. Nota anche che nell’implementazione non è necessario memorizzare tutta la successione x(k) basta solo operare su due vettori x0 e x1.


Io nn ho capito alcune cose, ad esempio: perchè nell'algoritmo dell'esercizio 1, c'è A(i,k) e non L ed U??
Come si possono fare le sommatorie in matlab??

grazie anticipatamente

[gandelf]

se qualcuno avesse una qualsiasi minima implementazione del secondo problema gliene sarei davvero grato..

[in_me_i_trust]

per jacobi detta A la matrice ,b il vettore termini noti,x_0 un qualunque vettore di partenza(non influenza la convergenza del metodo),itermax il numero massimo di iterazioni e toll la tolleranza dell'errore si potrebbe costruire cosi:
D=diag(diag(A));
D_1=inv(D);
C=A-D;
M=-D_1*C;
rho_j=max(abs(eig(M)));
if rho_j > 1
error('il raggio di convergenza è maggiore di 1');
end
for i=1:itermax
x=M*x_0+D_1*b;
if norm(x-x_0) <= toll*norm(x);
break
end
x_0=x;
end

per gauss seidel è praticamente uguale basta che ragioni un po con tril e triu..prova prima a scriverlo sotto forma di sommatorie o matrici dipende come te lo hanno spiegato..se prprio non ti riesce dimmelo

ciao!

[gandelf]

k grazie mille ora provo a ragiornarci su

[gandelf]

allora, quel che sono riuscito ad implementare è questo:

clc
clear

n=input('inserire n: ');
acc=input('inserire accuratezza: ');

A1=[1,-2,2;-1,1,-1;-2,-2,1];
A2=[2,-1,1;2,2,2;-1,-1,2];
B=(1/n)*rand(n);
DD=diag(diag(B));
B=B-DD;
A3=DD-DD*(B);

Ma da questo punto nn riesco piu ad andare avanti..

L'algoritmo che mi viene chiesto è questo:

For i:1,....,n
[x(i)]^(k+1) = ( b(i) - sum[per j=1 a i-1](A(ij) * x(j))^(k+1) - sum[per j=1+1 a n](A(ij) * x(j))^(k) ) 

Ho trovato delle soluzioni, ma nn si accostano manco lontanamente all'algoritmo qui sopra esposto..

ad esempio:

function [x; iter] = jacobi(a; b; x0; nmax; toll);
[n; n] = size(a); iter = 0;
r = b - a * x0; res = norm(r);
xold = x0; x = x0;%x vettore colonna
while res > toll & iter < nmax
iter = iter + 1;
for i = 1 : n
s = 0;
for j = 1 : i - 1; s = s + a(i; j) * xold(j); end
for j = i + 1 : n; s = s + a(i; j) * xold(j); end
x(i) = (b(i) - s)=a(i; i);
end
r = b - a * x; res = norm(r); xold = x;
end

oppure:

function [u,normz,k,D,L,U]=jacobi(A,b,u0,toll,nmax)
   %
   %
   % Metodo di Jacobi per la risoluzione di Au=b
   %
   % A=D-L-U
   
   % Parte diagonale
   if issparse(A) %per matrici sparse
    n=max(size(A));
    D=spdiags(spdiags(A,0),0,n,n);
   else %per matrici piene
    D=diag(diag(A));
   end
   % Parte triangolare inferiore (esclusa la diagonale principale)
   L=-tril(A,-1);
   % Parte triangolare superiore (esclusa la diagonale principale)
   U=-triu(A,1);
   
   
   % inizializzazione
   P=D;
   u=u0;
   
   % iterazioni di Jacobi
   for k=1:nmax
    r=b-A*u;
    z=P\r;
    u=u+z;
    normz(k)=norm(z,inf);
    if (normz(k)<=toll)
    break
    end
   end
   
   if (k==nmax)
    warning('Jacobi non e'' arrivato a convergenza!')
   else
    fprintf('\nJacobi e'' arrivato a convergenza in %d iterazioni\n',k);
   end
   
   return

come si può implementare quell'algoritmo? qualcuno ha qualche idea??

[gandelf]

la nuova versione dell'algoritmo è questa:

clc
clear

n=input('inserire n: ');
acc=input('inserire accuratezza: ');

A1=[1,-2,2;-1,1,-1;-2,-2,1];
A2=[2,-1,1;2,2,2;-1,-1,2];
B=(1/n)*rand(n);
DD=diag(diag(B));
B=B-DD;
A3=DD-DD*(B);

b=rand(3,1);

xold=zeros(size(b));

% costruzione matrice diagonale D
D = diag(diag(A)) ;

% costruzione DINV (inversa della matrice D)
DINV = inv(D) ;

% costruzione matrice B
B = eye(size(A)) - DINV*A ;

% costruzione vettore q
q = DINV * b ;

% inizializzazione del contatore delle iterazioni
iter = 0 ;

% inizializzazione dell'errore 
err  = 1 ;

% definizioni del numero massimo di iterazioni ammesse
maxiter = 100 ;

while (iter<=maxiter) & (err>acc)

% aggiornamento contatore delle iterazioni
   iter = iter + 1 ;
   
   % calcolo xnew (vettore x al passo iter)
   xnew = B * xold + q ;
   
   % calcolo dell'errore 
   err = norm(xnew-xold) ;
   
   % visualizzazione del numero dell'iterazione e dell'errore 
   display([ iter err ])
   
   % aggiornamento vettore xold
   xold = xnew ;
  
end

% visualizzazione del vettore soluzione calcolato
display('il vettore soluzione x è:')
x = xnew

% visualizzazione del vettore soluzione esatto
display('la soluzione esatta è:')
xesatto = A\b

% visualizzazione dell'errore commesso (differenza x-xesatto)
display('il vettore errore è:')
epsilon = abs(xesatto-x)
    

sn riuscito a collegare qualcosa ma continuo a nn capire cosa sia questo vettore q

[in_me_i_trust]

il tuo q deriva dal fatto che tu hai un sistema del tipo

D*x(k+1)=-(A-D)*x(k)+b

e quindi moltiplicando per DINV che è l'inversa di D hai :

x(k+1)=-DINV*(A-D)*x(k)+DINV*b

e detto DINV*b=q,si ha svolgendo i calcoli proprio l'espressione che scrivi tu (ma che coincide anche con quella che ti ho scritto io)

x(k+1)=(I-DINV*A)x(k)+q

cioè

x(k+1)=B*x(k)+q

per capire meglio come si arriva all'espressione da cui sono partito ti invito a guardare queste dispensine della mia prof; vai su

http://www.de.unifi.it/anum/morandi/morandi_didatt.htm

e clicca

comp.calc.num_05.pdf

nelle prime pagine trovi gauus seidel e jacobi + metodi di rilassamento ecc..ecc..

[gandelf]

in_me_I_trust devo assolutamente fare i complimenti alla tua prof e ringraziar te x avermi dato l'occasione di leggere queste slides.. sono fatte benissimo!

ora le leggo e poi magari se nn capisco qualcosa ti faccio sapere, grazie mille di nuovo

[gandelf]

ps ho trovato anche le slides 4.. x caso esistono anche le 1 2 e 3?

[in_me_i_trust]

non ce ne sono altre comunque nei complementi c'è tutto quello che ti serve

ciao^^