560 likes | 727 Views
5 a -6 a lezione di laboratorio. Laurea Ingegneria CIVILE Lauree Specialistiche in Ingegneria CHIMICA, ELETTRONICA, AMBIENTE. a.a. 2007-2008. Come creare un grafico 2-D. Sintassi per disegnare una curva con: specifica dei dati nel vettore x e/o y
E N D
5a-6a lezione di laboratorio Laurea Ingegneria CIVILE Lauree Specialistiche in Ingegneria CHIMICA, ELETTRONICA, AMBIENTE a.a. 2007-2008
Come creare un grafico 2-D Sintassi per disegnare una curva con: • specifica dei dati nel vettore x e/o y • specifica del colore e dello stile della linea plot(x, y, ‘colore_stilelinea’) Sintassi per disegnare più curve: plot(x1, y1,’r*’,x2, y2,‘b-’,...)
Completamento di un grafico Per completare un grafico si può aggiungere: • un titolo:title('Grafico ed …'); • la griglia:grid; • le label sugli assi: xlabel('x');ylabel('y')
Editore grafico Uno dei modi per migliorare l’aspetto di un grafico è il seguente: selezionare View e, dal menu che compare, scegliere Plot Edit Toolbar; nella finestra compare una seconda barra.
Le icone indicate dalle linee, consentono di inserire testo, frecce e linee su una figura. Cliccando sull’icona si ottengono i tool del plot che consentono di modificare la grafica della figura. Come operare con l’editor grafico
Come inserire un testo sulla figura • Selezionare l’icona • T dalla barra oppure • cliccare su ‘Insert’e • scegliereTextBox • Posizionarsi nel punto desiderato e scrivere • Scrivere x_3 per ottenere x3 (opp. x^3 perx3) • cliccare fuori dal riquadro per rendere attivo lo scritto
Come inserire una freccia sulla figura • Selezionare • l’icona • oppure dal menu di • Insert selezionare • Arrow • Posizionarsi nel • punto di inizio della • freccia, trascinare il • mouse tenendo • premuto il suo • tasto sinistro fino al • punto di arrivo della • freccia.
Esercizio 1 a - Scrivere un file script che consenta di disegnare, sull’intervallo [0,4], le funzioni: y=3*sin(pi*x) e y=exp(-0.2*x) nella stessa finestra grafica. Si consideri la partizione x=0:0.02:4. b - Inserire le label per gli assi x, y ed il titolo. c - Usare gtext per indicare i vari punti di intersezione delle due curve. d - Memorizzare il file col nome grafico.
File grafico.m clear all x=0:0.02:4; y=3*sin(pi*x); plot(x,y,'r'),xlabel('x');ylabel('y');grid hold on y1=exp(-0.2*x); plot(x,y1,'g') %osservare il numero delle intersezioni gtext('x1');gtext('x2');gtext('x3');gtext('x4'); % oppure %C=strvcat(‘x1’,’x2’,’x3’,’x4’); %gtext(C) title('Grafico ed intersezioni di 3*sin(pi*x) e exp(-0.2*x)') hold off
titolo etichetta asse y gtext etichetta asse x Risultato esercizio 1 » title('Grafico ed …') »gtext('x3') » ylabel('y') » xlabel('x')
Esercizio 2 Scrivere un file script che consenta di disegnare, nell’intervallo [-2,2] e su due “finestre grafiche distinte ”, il grafico della funzione: f(x)=exp(-x2)cos(20x) che viene definita nella function fun. Siutilizzino i due comandi MATLAB: plot per la figura 1 fplotper la figura 2. N.B.Nell’utilizzare il comando plot si può considerare la partizionex=[-2:0.1:2].
Soluzione esercizio 2: comando plot figure(1) x=[-2:0.1:2]; y=fun(x); plot(x,y),title('Comando plot') xlabel('x');ylabel('y'); grid functiony=fun(x) y=exp(-x.^2).*cos(20*x);
Soluzione esercizio 2: comando fplot figure(2) I=[-2,2]; fplot('fun',I), grid % Matlab ver.6 % utilizzabile anche nella ver. 7 title('Comando fplot') xlabel('x');ylabel('y') figure(3) % metodo alternativo I=[-2,2]; fplot(@fun,I), grid % Matlab ver.7 title('Comando fplot') xlabel('x');ylabel('y') functiony=fun(x) y=exp(-x.^2).*cos(20*x);
Comandi plot e fplot E’ possibile utilizzare i comandi plot e fplot, senza definire un file function esterno: x=[-2:0.1:2]; y=exp(-x.^2).*cos(20*x); figure(1) plot(x,y),grid %oppure f='exp(-x.^2).*cos(20*x)';%stringa %y=eval(f); %crea il vettore di dimensione % =length(x) %plot(x,y),grid f='exp(-x.^2).*cos(20*x)';% stringa figure(2) fplot(f,[-2,2]),grid
Comando ezplot Il comando ezplot consente di graficare funzioni date in forma implicita ed in forma parametrica figure(1) ezplot('x.^2-y.^2-1',[-4 4 -4 4]) %forma implicita
figure(2) ezplot('cos(t)','sin(t)') % forma parametrica
Più grafici in una finestra grafica E’ possibile inserire più di un grafico nella stessa finestra grafica . % Comando fplot fplot(@(x)[sin(2*x) x.*exp(-x)], [-1 1]) %Matlab7 1) % Comando plot x=-1:0.1:1; y1=sin(2*x); y2=x.*exp(-x); plot(x,y1,x,y2) 2) Alternativamente si può utilizzare il comando hold on: fplot('sin(2*x)',[-1 1]) hold on fplot('x.*exp(-x)',[-1 1]) hold off
Più finestre grafiche in una sola figura: subplot %file figure x=0.1:.1:5; subplot(2,3,1);plot(x,x); title('y=x');xlabel('x'); ylabel('y'); subplot(2,3,2);plot(x,x.^2); title('y= x^2');xlabel('x'); ylabel('y'); subplot(2,3,3),plot(x,x.^3); title('y= x^3');xlabel('x'); ylabel('y'); subplot(2,3,4),plot(x,cos(x)); title('y=cos(x)');xlabel('x');ylabel('y'); subplot(2,3,5),plot(x,cos(2*x)); title('y=cos(2x)');xlabel('x');ylabel('y') subplot(2,3,6),plot(x,cos(3*x)); title('y=cos(3x)');xlabel('x');ylabel('y') % I colori e lo spessore sono stati % aggiunti utilizzando i tool del plot
Comandi: plot e semilog Due modi di graficare la funzione Scala lineare su entrambi gli Assi x=0:0.1:50; y=exp(x); subplot(1,2,1);plot(x,y); grid;title('plot') subplot(1,2,2); semilogy(x,y); grid;title('semilogy') Scala logaritmica sull’Asse Y
Esercizio 3 (Esame 02/12/2002) Si considerino i sistemi lineari Ai xi=bi, i=1,2,3, con i vettori dei termini noti bi, i=1,2,3,scelti in modo che la soluzione dei sistemi sia i=[1,1,1,1]T, i=1,2,3. Supponiamo che:
Quesito 1 Si determini, mediante MATLAB, il condizionamento K2(Ai ), i = 1, 2, 3e si verifichi che K2(Ai) = (K2 ( A1 ))i , i=2, 3. Si spieghi il motivo di tale relazione e se ne prevedano le conseguenze.
Quesito 2a Si costruisca un file MATLAB: Cognome_Nome.m, che una volta avviato: a) faccia visualizzare una schermata con i dati personali (cognome, nome, matricola, corso di Laurea) ed una breve presentazione del problema;
Quesito 2b b) mediante un ciclo for, determini i dati Ai ,bi , i=1, 2, 3; risolva quindi i sistemi applicando il metodo di Gauss con pivoting parziale; calcoli l’errore relativo in norma 2;
Quesito 2c c) faccia visualizzare una tabella riassuntiva che riporti: intestazione: indice iter soluzione errore e su ogni riga il valore dell’indice dellamatrice, il numero di iterazioni effettuate nel raffinamento, la soluzione corrispondente scritta come vettore riga e l’errore relativo ; Si utilizzino i seguenti formati di stampa: 1 cifra intera per il valore di ; 2 cifre intere per il valore di ; 10 cifre decimali e formato virgola fissa per le componenti di ; 2 cifre decimali e formato esponenziale per .
Soluzione teorica del Quesito 1 Proprietà della matrice e conseguenze:
K2(A1) Soluzione teorica del Quesito 1 Proprietà delle matrici per la simmetria di quindi: Analogamente per i = 3.
Istruzioni relative al Quesito 1 % file script: punto1.m clear all disp('Numero di condizionamento delle matrici Ai') A1=[15 6 8 11 ; 6 6 5 3 ; 8 5 7 6; 11 3 6 9]; cond_Ai=[];cond_A=[]; for i =1:3 Ai=A1^i; cond_Ai=[cond_Ai,cond(Ai)]; % vettore dei cond(Ai) cond_A=[cond_A,cond(A1)^i]; % vettore dei cond(A1)^i end disp('cond(Ai)') disp(num2str(cond_Ai,'%13.3e')) disp('(cond(A1))^i') disp(num2str(cond_A,'%13.3e'))
Output punto1 >> punto1 Numero di condizionamento delle matrici Ai cond(Ai) 6.499e+003 4.224e+007 2.746e+011 (cond(A1))^i 6.499e+003 4.224e+007 2.745e+011 Indicheremo con K1, K2, K3 il condizionamento in norma 2 delle matrici A1, A2, A3 rispettivamente. Conseguenze del numero di condizionamento grande?
Numero di cifre significative perse >> nc=log10(K2) nc = 7.6257 10 cifre signific., 9 decimali corretti Calcolo della soluzione di >> x2=A2\b2 % Operatore \ x2 = 1.00000000037954 0.99999999968201 1.00000000023455 0.99999999948547 >> err2=norm(x2-alpha)/norm(alpha) err2 = 3.76e-010=.376e-009
Oltre che sul numero di cifre significative che si perdono, come incide il valore grande di K2(A2)? >> A2m=A2; >> A2m(2,2)=A2(2,2)+1e-3; % perturbazione data % sulla matrice >>pert=norm(A2-A2m)/norm(A2) pert = 1.0546e-006 % entità della perturbazione >> x2m=A2m\b2 % soluzione perturbata x2m = 2.06050813471394 0.11143302866360 1.65544907516937 -0.43770893706830 >> err2m=norm(x2m-alpha)/norm(alpha) err2m = 1.0501>100% !!!
Istruzioni relative al Quesito 2a clear all disp('Cognome e nome studente: XXXX XXX') disp('Numero di matricola: XXXX') disp('Corso di Laurea: XXXX') disp(' ') disp('Questo programma calcola e visualizza la soluzione dei ') disp('sistemi lineari A_i x_i=b_i, i =1,2,3, con i vettori b_i tali che sia') disp('alpha=[1,1,1,1]'',essendo: ') A1=[15 6 8 11;6 6 5 3;8 5 7 6;11 3 6 9]; disp('A1=');disp(A1) disp( 'e A_i= A1^i per i=2,3.')
Output file Cognome_Nome.m >> Cognome_Nome Cognome e nome studente: XXXX XXX Numero di matricola: XXXX Corso di Laurea:XXXX Questo programma calcola e visualizza la soluzione dei sistemi lineari A_i x_i=b_i, i =1,2,3, con i vettori b_i tali che sia alpha=[1 1 1 1]', essendo: A1= 15 6 8 11 6 6 5 3 8 5 7 6 11 3 6 9 e A_i= A1^i per i=2,3.
Istruzioni di Gausspv_r.m Istruzioni relative al Quesito 2b tab=[]; toll=1e-13; alpha=ones(4,1); % soluzione for i =1:3 A=A1^i; b=A*alpha; % vettore termini noti % [L,U,P] = lu(A); %y=L\(P*b); %x=U\y; %[x,iter]=Raff_iter(A,b,L,U,P,x,toll); residuo=b-A*x; norm_residuo=norm(b-A*x); err=norm(alpha-x)/norm(alpha); tab=[tab;[i,iter,x',err,norm_residuo]]; end %La chiamata di Gausspv_r è: [x,iter] = Gausspv_r(A,b,toll)
Numero di cifre significative perse >> nc=log10(K3) nc = 11.4386 Nell’ultimo caso calcolando la stima del numero di cifre che si perdono, si ottiene: Risultati File Quesito 2b Tabella Quesito 2c fprintf('i iter \t\t\t soluzione \t\t \t\t errore residuo\n') fprintf('%1d %2d %14.10f %14.10f %14.10f %14.10f %10.2e %10.2e \n',tab'); i iter soluzione errore residuo 1 0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 8.17e-015 3.55e-015 2 0 1.0000000007 0.9999999994 1.0000000004 0.9999999991 6.81e-010 1.97e-013 3 0 0.9999998145 1.0000001554 0.9999998854 1.0000002514 1.84e-007 0.00e+000
Istruzioni utilizzate in Raff_iter function [x,iter]=Raff_iter(A,b,L,U,P,x,toll) %_________________________________________ ... iter=0; residuo=b-A*x; while norm(residuo)>toll*norm(b)& iter<100 y=L\(P*residuo); err=U\y; x=x+err; residuo=b-A*x; iter=iter+1; end if iter==100 disp('Raggiunto il numero massimo di iterazioni') end
Esercizio 4 Sia dato il sistema lineare avente la matrice dei coefficienti ed il vettore dei termini noti così assegnati:
Quesiti 1 e 2 • Si studi la convergenza dei metodi di Jacobi, Gauss-Seidel e Rilassamento in serie (SOR) per il sistema assegnato. • Si dica quale di questi metodi è il più veloce giustificando teoricamente la risposta e calcolando, mediante Matlab, i raggi spettrali delle rispettive matrici di iterazione ed il valore ottimale del parametro per il metodo SOR.
Quesito 3a 3. Si costruisca un file MATLAB che: • calcoli la soluzione numerica del problema assegnato applicando il metodo con convergenza migliore e fissando una precisione non inferiore a 1.e-4, nmax=20 ed un vettore di innesco pari a x0=[1 1 1 1 1 1]T;
Quesito 3b b) faccia visualizzare una tabella riassuntiva che riporti: intestazione: iterazioni soluzione residuo; e su ogni riga, il numero dell’ iterazione , la soluzione approssimata corrispondente e la norma del residuo ; Si utilizzino i seguenti formati di stampa: 2 cifre intere per il valore di ; 5 cifre decimali e formato virgola fissa per la soluzione approssimata ; 1 cifra decimale e formato esponenziale per la norma del residuo .
Convergenza dei metodi • Caratteristiche di A: • diagonalmente dominante Jacobi conv. • tridiagonale anche Gauss-Seidel converge; inoltre: c) simmetrica d) definita positiva perché è a), c) e tutti gli elementi sulla diagonale principale sono positivi.
Qual è il metodo più veloce? Per la proprietà d) SOR converge per è inoltre il metodo più veloce se si assume:
Istruzioni relative al Quesito 2 I0=[4 5 6 5 3 2];I1=[1 2 3 -1 1]; A=diag(I0)+diag(I1,-1)+diag(I1,1) D=diag(diag(A)); [n,m]=size(A); B_J=eye(n)-inv(D)*A; % metodo di Jacobi rho_J=max(abs(eig(B_J))) R_J=-log(rho_J) omega=1; % metodo di Gauss-Seidel OE=omega*tril(A,-1); B_GS=eye(n)-omega*inv(D+OE)*A; rho_GS=max(abs(eig(B_GS))) R_GS=-log(rho_GS) omega_ott=2/(1+sqrt(1-rho_GS)) % metodo SOR OE=omega_ott*tril(A,-1); B_r=eye(n)-omega_ott*inv(D+OE)*A; rho_r=max(abs(eig(B_r))) R_r=-log(rho_r)
Risultati file Quesito 2 rho_J = 0.7196R_J = 0.3290 rho_GS = 0.5179R_GS = 0.6580 omega_ott = 1.1804 rho_r = 0.1804R_r = 1.7126
Istruzioni relative al Quesito 3 b=[3 -2 1 -1 3 4]';K=cond(A,inf); precisione=input('precisione = ');% 1.e-4 toll=precisione/K% toll = 9.3310e-006 x0=ones(n,1);omega=omega_ott;nmax=20; [x,iter,res,rho]=Gauss_Seidel_ril(A,b,x0,omega,nmax,toll); it=[0:iter]';tab=[it x res]; s='--------------------------------------------'; disp(s) fprintf('iter soluzione errore\n') fprintf('%2d %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %9.1e\n',tab');
Risultati file Quesito 3 • iter soluzione errore • ------------------------------------------------------------------ • 0 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 2.5e+000 • 1 0.40980 -1.22145 -0.09326 -0.11434 0.56155 1.84897 1.1e+000 • 2 1.17181 -0.48443 0.47164 -0.41692 0.18755 1.91655 4.2e-001 • 3 0.81686 -0.80030 0.67260 -0.59295 0.15916 1.92112 1.3e-001 • 4 0.97411 -0.87533 0.76977 -0.63672 0.14526 1.92850 3.4e-002 • 5 0.96788 -0.90620 0.79022 -0.64659 0.14099 1.92969 8.7e-003 • 6 0.97811 -0.91270 0.79491 -0.64914 0.14028 1.92989 1.1e-003 • 7 0.97819 -0.91376 0.79599 -0.64961 0.14015 1.92993 3.4e-004 • 8 0.97849 -0.91415 0.79622 -0.64972 0.14011 1.92995 5.1e-005 • 9 0.97855 -0.91421 0.79627 -0.64974 0.14010 1.92995 1.3e-005 • 0.97855 -0.91422 0.79628 -0.64975 0.14010 1.92995 2.5e-006 Si invitano gli studenti a determinare la soluzione con precisione = 1.e-8 e riportare le componenti della soluzione approssimata con 10 cifre e virgola fissa.
Esercizio 5 (Esame 05-12-05) Sia data la seguente matrice: 0 evettore b tale che la soluzione del sistemasia
Quesiti 1, 2 e 3 1 - Dopo aver determinato con MATLAB gli autovalori, si deduca motivando la risposta, la caratteristica fondamen_ tale della matrice; 2 - si determini il condizionamento in norma 2 e si dica, sempre motivando la risposta, se il sistema è ben condi_ zionato calcolando, inoltre, il numero di cifre significative che si perdono, rispetto alle 16 del MATLAB, risolvendo il sistema; 3 - Si costruisca un file MATLAB: Cognome_NomeStudente.m che, una volta avviato: