Introduzione al MatLab

(7)

1

Metodi Numerici con Elementi di Programmazione

A.A. 2019-2020

2

Docente: Vittoria Bruni Email: vittoria.bruni@sbai.uniroma1.it Ufficio: Via A. Scarpa, Pal. B, I piano, Stanza n. 16 Tel. 06 49766648 Ricevimento: consultare la sezione Avvisi nella pagina web dedicata al corso

3

Il materiale didattico è disponibile sul sito http://ingaero.uniroma1.it/

nella pagina dedicata al corso Metodi Numerici con Elementi di

Programmazione

Per consultazione: Getting Started with MatLab – The mathworks

www.mathworks.com

Le oscillazioni di un pendolo semplice possono essere descritte dall'equazione differenziale del secondo ordine non lineare dove L è la lunghezza del pendolo, g è l'accelerazione di gravità e è l'angolo tra il pendolo e la verticale. Il problema è completato con le condizioni iniziali Scrivere tre funzioni Matlab per la soluzione di sistemi di equazioni differenziali del primo ordine rispettivamente mediante il metodo di Eulero, di Heun e di Runge Kutta del IV ordine con passo di discretizzazione h. Usare le funzioni per approssimare la soluzione del problema del pendolo semplice non lineare nei primi 10 secondi di oscillazione sapendo che Confrontare i risultati ottenuti con quelli prodotti dalla funzione Matlab ode45.m Confrontare i risultati al variare di L, del tempo di osservazione e delle condizioni iniziali.

Esercizio

Soluzione L'equazione differenziale del secondo ordine non lineare si trasforma nel seguente sistema di equazioni differenziali del primo ordine

1. Si scrive la funzione Matlab che descrive il problema specifico

function [f] = eqPendolo(t,theta) % Equazione differenziale che descrive le oscillazioni del pendolo % SINTASSI % f = eqPendolo(t,theta) % % INPUT % t : Valore della variabile indipendente % theta : Valore della funzione vettoriale theta % % OUTPUT %f: Valore di f(t, theta) in t e theta % costante di gravità g = 9.81; % lunghezza del pendolo L = 0.6; % equazione differenziale f = [theta(2);-g/L*sin(theta(1))];

2. Si scrivono le funzioni Matlab che risolvono un sistema di equazioni differenziali rispettivamente con i metodi di Eulero, Heun e Runge Kutta del IV ordine

function [T] = Eulero_sistemi(I,Y0,h,ffun) % Soluzione numerica del sistema di n equazioni differenziali del primo ordine % Y’=ffun(x,Y) prodotta mediante il metodo di Eulero esplicito % INPUT % I = intervallo in cui si cerca la soluzione [x0, xN] % Y0= vettore delle condizioni iniziali % h = passo di discretizzazione % ffun = funzione del problema differenziale % % OUTPUT % T = matrice contenente la soluzione [xi;Yi+ nell’intervallo *x0,xN]. % La prima riga di T contiene i nodi in cui stata stimata la soluzione % Le righe successive contengono le stime della funzione soluzione del % problema differenziale e delle sue derivate fino all’odine n-1 n_step = (I(end)-I(1))/h +1; n=length(Y0); xi(1) = I(1); Yi(1:n,1) = Y0(:);

for i = 2:n_step x = xi(i-1); Y = Yi(:,i-1); xi(i) = x + h; Yi(:,i) = Y + h*ffun(x,Y); end T = [xi;Yi];

function [T] = Heun_sistemi(I,Y0,h,ffun) % Soluzione numerica del sistema di n equazioni differenziali del primo ordine % Y’=ffun(x,Y) prodotta mediante il metodo di Heun (secondo ordine) % INPUT % I = intervallo in cui si cerca la soluzione [x0, xN] % Y0= vettore delle condizioni iniziali % h = passo di discretizzazione % ffun = funzione del problema differenziale % % OUTPUT % T = matrice contenente la soluzione [xi;Yi+ nell’intervallo *x0,xN+. % La prima riga di T contiene i nodi in cui stata stimata la soluzione % Le righe successive contengono le stime della funzione soluzione del % problema differenziale e delle sue derivate fino all’odine n-1 n_step = (I(end)-I(1))/h +1; n = length(Y0); xi(1) = I(1); Yi(1:n,1) = Y0(:);

for i = 2:n_step x = xi(i-1); Y = Yi(:,i-1); K1 = ffun(x,Y); K2 = ffun(x+h,Y+h*K1); xi(i) = x + h; Yi(:,i) = Y + h*(K1 + K2)/2; End T = [xi;Yi];

function [T] = RK_sistemi(I,Y0,h,ffun) % Soluzione numerica del sistema di n equazioni differenziali del primo ordine % Y’=ffun(x,Y) prodotta mediante il metodo di Runge Kutta del quarto ordine % INPUT % I = intervallo in cui si cerca la soluzione [x0, xN] % Y0= vettore delle condizioni iniziali % h = passo di discretizzazione % ffun = funzione del problema differenziale % % OUTPUT % T = matrice contenente la soluzione [xi;Yi+ nell’intervallo *x0,xN+. % La prima riga di T contiene i nodi in cui stata stimata la soluzione % Le righe successive contengono le stime della funzione soluzione del % problema differenziale e delle sue derivate fino all’odine n-1 n_step = (I(end)-I(1))/h +1; n=length(Y0); xi(1) = I(1); Yi(1:n,1) = Y0(:);

for i = 2:n_step x = xi(i-1); Y = Yi(:,i-1); K1 = ffun(x,Y); K2 = ffun(x+0.5*h,Y+0.5*h*K1); K3 = ffun(x+0.5*h,Y+0.5*h*K2); K4 = ffun(x+h,Y+h*K3); xi(i) = x + h; Yi(:,i) = Y + h*(K1+2*K2+2*K3+K4)/6; end T = [xi;Yi];

3. Si scrive uno script Matlab che risolve richiama le funzioni Eulero_sistemi , Heun_sistemi e RK_sistemi passando gli input relativi al problema differenziale specifico

% File mainPendolo.m close all clear all I = [0, 10]; Y0 = [pi/6;0]; h = 0.05; [T_E] = Eulero_sistemi(I,Y0, h,@eqPendolo); [T_H] = Heun_sistemi(I,Y0, h,@eqPendolo); [T_RK] = RK_sistemi(I,Y0, h,@eqPendolo);

% grafico dei risultati del metodo di Eulero figure(1) subplot(2,1,1) plot(T_E(1,:),T_E(2,:),’o’) xlabel(‘t’) ylabel(‘\theta(t) - Eulero’) title(‘oscillazione’) subplot(2,1,2) plot(T_E(1,:),T_E(3,:),’*’) xlabel(‘t’) ylabel(‘d\theta(t) - Eulero’) title(‘velocita’’ angolare’) pause

4. Si visualizzano i risultati ottenuti con i tre metodi in tre finestre grafiche distinte: si visualizzano separatamente l’angolo dell’oscillazione (posizione del pendolo rispetto alla verticale) e la velocità angolare (derivata dell’angolo di oscillazione)

% grafico dei risultati del metodo di Heun figure(2) subplot(2,1,1) plot(T_H(1,:),T_H(2,:),’o’) xlabel(‘t’) ylabel(‘\theta(t) - Heun’) title(‘oscillazione’) subplot(2,1,2) plot(T_H(1,:),T_H(3,:),’*’) xlabel(‘t’) ylabel(‘d\theta(t) - Heun’) title(‘velocita’’ angolare’) pause

% grafico dei risultati del metodo di Runge-Kutta IV ordine figure(3) subplot(2,1,1) plot(T_RK (1,:),T_RK (2,:),’o’) xlabel(‘t’) ylabel(‘\theta(t) – RK4’) title(‘oscillazione’) subplot(2,1,2) plot(T_RK (1,:),T_RK (3,:),’*’) xlabel(‘t’) ylabel(‘d\theta(t) – RK4’) title(‘velocita’’ angolare’) pause

% grafico dei risultati del metodo di Eulero figure(4) for i=1:length(T_E(1,:)) plot(T_E(2,1:i),T_E(3,1:i)) xlabel('rad'), ylabel('rad/s') axis([min(T_E(2,:)) max(T_E(2,:)) min(T_E(3,:)) max(T_E(3,:))]) title(*‘Eulero -- t=' num2str(T_E(1,i))]) hold on plot(T_E(2,i),T_E(3,i),'o','markersize',10,'markerFaceColor',‘r','markerEdgeColor',‘r') pause(.01) frame = getframe(4); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,'E_simple.gif','gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,'E_simple.gif','gif','DelayTime',0,'writemode','append'); end hold off end, pause

4b. Per ognuno dei tre metodi, si visualizza il grafico (angolo di oscillazione, velocità angolare) con animazione del risultato

% grafico dei risultati del metodo di Heun figure(5) for i=1:length(T_H(1,:)) plot(T_H(2,1:i),T_H(3,1:i)) xlabel('rad'), ylabel('rad/s') axis([min(T_H(2,:)) max(T_H(2,:)) min(T_H(3,:)) max(T_H(3,:))]) title(*‘Heun -- t=' num2str(T_H(1,i))]) hold on plot(T_H(2,i),T_H(3,i),'o','markersize',10,'markerFaceColor',‘r','markerEdgeColor',‘r') pause(.01) frame = getframe(5); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,‘H_simple.gif','gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,‘H_simple.gif','gif','DelayTime',0,'writemode','append'); end hold off end, pause

4b. Per ognuno dei tre metodi, si visualizza il grafico (angolo di oscillazione, velocità angolare) con animazione del risultato

% grafico dei risultati del metodo di Runge-Kutta IV ordine figure(6) for i=1:length(T_RK(1,:)) plot(T_RK (2,1:i),T_RK (3,1:i)) xlabel('rad'), ylabel('rad/s') axis([min(T_RK(2,:)) max(T_RK(2,:)) min(T_RK(3,:)) max(T_RK(3,:))]) title(*‘RK4 -- t=' num2str(T_RK (1,i))]) hold on plot(T_RK (2,i),T_RK(3,i),'o','markersize',10,'markerFaceColor',‘r','markerEdgeColor',‘r') pause(.01) frame = getframe(6); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,‘RK_simple.gif','gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,‘RK_simple.gif','gif','DelayTime',0,'writemode','append'); end hold off end, pause

4b. Per ognuno dei tre metodi, si visualizza il grafico (angolo di oscillazione, velocità angolare) con animazione del risultato

5. Usare lo help di Matlab per conoscere cosa implementa la funzione ode45 e come usarla;

% File mainPendolo.m options = odeset('Vectorized','on','InitialStep',h); [t,Y]=ode45(@eqPendolo,I,Y0,options);

La funzione Matlab ode45 implementa una coppia

di metodi Runge-Kutta espliciti di ordine 4 e 5 per definire un passo

di discretizzazione adattivo.

L’algoritmo usato implementa il metodo di Dormand–Prince (RKDP) e

il passo è scelto in modo da migliorare la stima del metodo del

quarto ordine

Usare help!!!

5. Usare lo help di Matlab per conoscere cosa implementa la funzione ode45 e come usarla; quindi, completare lo script mainPendolo.m utilizzando la funzione ode45 per il problema specifico

% File mainPendolo.m options = odeset('Vectorized','on','InitialStep',h); [t,Y]=ode45(@eqPendolo,I,Y0,options); figure(7) subplot(2,1,1) plot(t,Y(:,1),’o’) xlabel(‘t’) ylabel(‘\theta(t) – ode45’) title(‘oscillazione’) subplot(2,1,2) plot(t,Y(:,2),’*’) xlabel(‘t’) ylabel(‘d\theta(t) – ode45’) title(‘velocita’’ angolare’) pause

% grafico dei risultati ode45 figure(8) for i=1:length(t) plot(Y(1:i,1),Y (1:i,2)) xlabel('rad'), ylabel('rad/s') axis([min(Y(:,1)) max(Y(:,1)) min(Y(:,2)) max(Y(:,2))]) title(*‘ode45 -- t=' num2str(t(i))]) hold on plot(Y (I,1),Y(I,2),'o','markersize',10,'markerFaceColor',‘r','markerEdgeColor',‘r') pause(.01) frame = getframe(8); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,‘ode45_simple.gif','gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,‘ode45_simple.gif','gif','DelayTime',0,'writemode','append'); end hold off end pause

5. Usare lo help di Matlab per conoscere cosa implementa la funzione ode45 e come usarla; quindi, completare lo script mainPendolo.m utilizzando la funzione ode45 per il problema specifico e confrontare il risultato con quello ottenuto con i metodi di Eulero, Heun e RK4

figure(3) subplot(2,1,1), hold on plot(t,Y(:,1),’ro’);, legend(‘Rksistemi’,’ode45’) subplot(2,1,2), hold on plot(t,Y(:,2),’r*’); legend(‘Rksistemi’,’ode45’) pause

Scrivere le istruzioni necessarie per confrontare il risultato di ode45 anche con il risultato prodotto dai metodi di Eulero e di Heun

6. Richiamare il file mainPendolo.m dal Command Window; osservare e commentare i risultati

Caso h = 0.05 tf=10 Eulero % dal Command window >> mainPendolo

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

Caso h = 0.05 tf=10 Heun

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

Caso h = 0.05 tf=10 RK4

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

Caso h = 0.05 tf=10

Caso h = 0.05 tf=10

7. Modificare i parametri del problema e/o dei metodi numerici, richiamare il file mainPendolo.m dal Command Window; osservare e commentare i risultati modificare valori dei parametri in File mainPendolo.m I = [0, 1]; Y0 = [pi/6;0]; h = 0.05; e richiamare lo script >> mainPendolo

EULERO

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

HEUN

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

RK4

Nota: Osservare la griglia di discretizzazione usata dalla funzione ode45

7. Modificare i parametri del problema e/o dei metodi numerici, richiamare il file mainPendolo.m dal Command Window; osservare e commentare i risultati

Caso h = 0.5 Eulero

Caso h = 0.5 Heun

Caso h = 0.5 RK4

7. Modificare i parametri del problema e/o dei metodi numerici, richiamare il file mainPendolo.m dal Command Window; osservare e commentare i risultati

Caso h = 0.1 L=0.06 tf=1 Eulero

Nota: In questo caso è necessario modificare sia valori dei parametri in mainPendolo.m che nella funzione eqPendolo.m

Caso h = 0.1 L=0.06 tf=1 Heun

Caso h = 0.1 L=0.06 tf=1 RK4

Esercizio

38

Per approssimare la soluzione di questo problema differenziale è necessario: - modificare la funzione eqPendolo in quanto è cambiata la funzione del problema differenziale e i parametri da cui dipende - richiamare lo script mainPendolo.m ne

1. Si scrive la funzione Matlab che descrive il problema specifico

function [f] = eqPendolo(t,theta) % Equazione differenziale che descrive le oscillazioni del pendolo % SINTASSI % f = eqPendolo(t,theta) % % INPUT % t : Valore della variabile indipendente % theta : Valore della funzione vettoriale theta % % OUTPUT %f: Valore di f(t, theta) in t e theta % costante di gravità g = 9.81; % lunghezza del pendolo L = 0.6; % coefficiente di smorzamento b=0.6; % b=k*L/m, con m= massa del pendolo, k=coefficiente di resistenza % equazione differenziale f = [theta(2);-g/L*sin(theta(1))-b*theta(2)];

1. Si scrive la funzione Matlab che descrive il problema specifico

function [f] = eqPendolo(t,theta) % Equazione differenziale che descrive le oscillazioni del pendolo % SINTASSI % f = eqPendolo(t,theta) % % INPUT % t : Valore della variabile indipendente % theta : Valore della funzione vettoriale theta % % OUTPUT %f: Valore di f(t, theta) in t e theta % costante di gravità g = 9.81; % lunghezza del pendolo L = 0.6; % coefficiente di smorzamento b=0.6;% b=k*L/m, con m= massa del pendolo, k=coefficiente di resistenza % equazione differenziale f = [theta(2);-g/L*sin(theta(1))-b*theta(2)];

Modificare in modo che i parametri siano variabili di

input

41

% File mainPendolo.m I = [0, 10]; Y0 = [pi/6;0]; h = 0.05; [T_E] = Eulero_sistemi(I,Y0, h,@eqPendolo); [T_H] = Heun_sistemi(I,Y0, h,@eqPendolo); [T_RK] = RK_sistemi(I,Y0, h,@eqPendolo);

2. Definire i valori iniziali nello script mainPendolo

2. Richiamare lo script mainPendolo dal Command Window

>> mainPendolo

42

43

44

45

Esercizio

1. Si scrive la funzione Matlab che descrive il problema specifico

function [f] = eqSpacecraft(t,V) % Equazione differenziale che descrive il moto della navicella spaziale lanciata ad % un'altitudine H su livello del mare ad una velocita v0 % SINTASSI % f = eqSpacecraft(t,V) % % INPUT % t : Valore della variabile indipendente % V: Valore della funzione vettoriale contenente le coordinate polari e le loro % derivate rispetto al tempo *r, theta, r’, theta’++ % % OUTPUT %f: Valore di f(t, V) in t e V % costante gravitazionale G = 6.672*10^-11; % massa della terra MT = 5.9742*10^24; % equazione differenziale f = [V(3); V(4); V(1)*V(4)^2-G*MT/(V(1)^2); -2*V(3)*V(4)/(V(1))];

2. Si scrive lo script in cui si definiscono i parametri, si chiamano le funzioni che implementano i metodi numerici, si confrontano i risultati, etc.

% File mainSpacecraft.m close all clear all nameE='Espace.gif'; nameH='Hspace.gif'; nameRK4='RK4space.gif'; nameode45='ode45space.gif'; namefun=@eqSpacecraft; I = [0, 1200]; H = 772*10^3; RT = 6378.14*10^3; v0 = 6700; r0 = RT+H; Y0 = [r0;0;0;v0/r0]; h = 50; [T_E] = Eulero_sistemi(I,Y0,h,namefun); [T_H] = Heun_sistemi(I,Y0,h,namefun); [T_RK] = RK_sistemi(I,Y0,h,namefun);

% grafico dei risultati del metodo di Eulero figure(1) subplot(2,2,1) plot(T_E(1,:),T_E(2,:),'o') xlabel('t') ylabel('r(t) - Eulero') title('r') subplot(2,2,2) plot(T_E(1,:),T_E(3,:),'*') xlabel('t') ylabel('theta(t) - Eulero') title('\theta') subplot(2,2,3) plot(T_E(1,:),T_E(4,:),'o') xlabel('t') ylabel('dr(t) - Eulero') title('dr(t)') subplot(2,2,4) plot(T_E(1,:),T_E(5,:),'*') xlabel('t') ylabel('d\theta(t) - Eulero') title('d\theta(t)') pause

% grafico dei risultati del metodo di Heun figure(2) subplot(2,2,1) plot(T_H(1,:),T_H(2,:),'o') xlabel('t') ylabel('r(t) - Heun') title('r') subplot(2,2,2) plot(T_H(1,:),T_H(3,:),'*') xlabel('t') ylabel('theta(t) - Heun') title('\theta') subplot(2,2,3) plot(T_H(1,:),T_H(4,:),'o') xlabel('t') ylabel('dr(t) - Heun') title('dr(t)') subplot(2,2,4) plot(T_H(1,:),T_H(5,:),'*') xlabel('t') ylabel('d\theta(t) - Heun') title('d\theta(t)') pause

% grafico dei risultati del metodo di Runge-Kutta IV ordine figure(3) subplot(2,2,1) plot(T_RK(1,:),T_RK(2,:),'o') xlabel('t') ylabel('r(t) - RK4') title('r') subplot(2,2,2) plot(T_RK(1,:),T_RK(3,:),'*') xlabel('t') ylabel('theta(t) - RK4') title('\theta') subplot(2,2,3) plot(T_RK(1,:),T_RK(4,:),'o') xlabel('t') ylabel('dr(t) - RK4') title('dr(t)') subplot(2,2,4) plot(T_RK(1,:),T_RK(5,:),'*') xlabel('t') ylabel('d\theta(t) - RK4') title('d\theta(t)') pause

% grafico dei risultati del metodo di Eulero figure(4) for i=1:length(T_E(1,:)) plot(T_E(2,1:i).*cos(T_E(3,1:i)),T_E(2,1:i).*sin(T_E(3,1:i))) xlabel('rcos(\theta)') ylabel('rsin(\theta)') axis([min(T_E(2,:).*cos(T_E(3,:))) max(T_E(2,:).*cos(T_E(3,:))) min(T_E(2,:).*sin(T_E(3,:))) max(T_E(2,:).*sin(T_E(3,:)))])

title(['Eulero -- t=' num2str(T_E(1,i))]) hold on plot(RT*cos(T_E(3,:)),RT*sin(T_E(3,:)),'k:') plot(T_E(2,i)*cos(T_E(3,i)),T_E(2,i)*sin(T_E(3,i)),'o','markersize',10,'markerFaceColor','r','markerEdgeColor','r') frame = getframe(4); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,nameE,'gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,nameE,'gif','DelayTime',0,'writemode','append'); end pause(.1) hold off end pause

% grafico dei risultati del metodo di Heun figure(5) for i=1:length(T_H(1,:)) plot(T_H(2,1:i).*cos(T_H(3,1:i)),T_H(2,1:i).*sin(T_H(3,1:i))) xlabel('rcos(\theta)') ylabel('rsin(\theta)') axis([min(T_H(2,:).*cos(T_H(3,:))) max(T_H(2,:).*cos(T_H(3,:))) min(T_H(2,:).*sin(T_H(3,:))) max(T_H(2,:).*sin(T_H(3,:)))])

title(['Heun -- t=' num2str(T_H(1,i))]) hold on plot(RT*cos(T_H(3,:)),RT*sin(T_H(3,:)),'k:') plot(T_H(2,i)*cos(T_H(3,i)),T_H(2,i)*sin(T_H(3,i)),'o','markersize',10,'markerFaceColor','r','markerEdgeColor','r')

pause(.1) frame = getframe(5); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,nameH,'gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,nameH,'gif','DelayTime',0,'writemode','append'); end hold off end pause

% grafico dei risultati del metodo di Runge-Kutta IV ordine figure(6) for i=1:length(T_RK(1,:)) plot(T_RK(2,1:i).*cos(T_RK(3,1:i)),T_RK(2,1:i).*sin(T_RK(3,1:i))) xlabel('rcos(\theta)') ylabel('rsin(\theta)') axis([min(T_RK(2,:).*cos(T_RK(3,:))) max(T_RK(2,:).*cos(T_RK(3,:))) min(T_RK(2,:).*sin(T_RK(3,:))) max(T_RK(2,:).*sin(T_RK(3,:)))])

title(['RK4 -- t=' num2str(T_RK (1,i))]) hold on plot(RT*cos(T_RK(3,:)),RT*sin(T_RK(3,:)),'k:') plot(T_RK(2,i)*cos(T_RK(3,i)),T_RK(2,i)*sin(T_RK(3,i)),'o','markersize',10,'markerFaceColor','r','markerEdgeColor','r')

pause(.1) frame = getframe(6); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,nameRK4,'gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,nameRK4,'gif','DelayTime',0,'writemode','append'); end hold off end pause

options = odeset('Vectorized','on','InitialStep',h); [t,Y]=ode45(namefun,I,Y0,options); figure(7) subplot(2,2,1) plot(t,Y(:,1),'o') xlabel('t') ylabel('r(t) – ode45') title('r') subplot(2,2,2) plot(t,Y(:,2),'*') xlabel('t') ylabel('\theta(t) – ode45') title('\theta(t)') subplot(2,2,3) plot(t,Y(:,3),'o') xlabel('t') ylabel('dr(t) – ode45') title('dr(t)') subplot(2,2,4) plot(t,Y(:,4),'*') xlabel('t') ylabel('d\theta(t) – ode45') title('d\theta(t)') pause

% grafico dei risultati ode45 figure(8) for i=1:length(t) plot(Y(1:i,1).*cos(Y(1:i,2)),Y(1:i,1).*sin(Y(1:i,2))) xlabel('rcos(\theta)') ylabel('rsin(\theta)') axis([min(Y(:,1).*cos(Y(:,2))) max(Y(:,1).*cos(Y(:,2))) min(Y(:,1).*sin(Y(:,2))) max(Y(:,1).*sin(Y(:,2)))])

title(['ode45 -- t=' num2str(t(i))]) hold on plot(RT*cos(Y(:,2)),RT*sin(Y(:,2)),'k:') plot(Y(i,1)*cos(Y(i,2)),Y(i,1)*sin(Y(i,2)),'o','markersize',10,'markerFaceColor','r','markerEdgeColor','r')

pause(.1) frame = getframe(8); img = frame2im(frame); [imgind,cm] = rgb2ind(img,256); if i == 1 imwrite(imgind,cm,nameode45,'gif','DelayTime',0,'loopcount',inf); else imwrite(imgind,cm,nameode45,'gif','DelayTime',0,'writemode','append'); end hold off end pause

% confronti con ode45 figure(1) subplot(2,2,1), hold on plot(t,Y(:,1),'ro');, legend('Eulerosistemi','ode45') subplot(2,2,2), hold on plot(t,Y(:,2),'r*'); legend('Eulerosistemi','ode45') subplot(2,2,3), hold on plot(t,Y(:,3),'r*'); legend('Eulerosistemi','ode45') subplot(2,2,4), hold on plot(t,Y(:,4),'r*'); legend('Eulerosistemi','ode45') pause

% confronti con ode45 figure(2) subplot(2,2,1), hold on plot(t,Y(:,1),'ro');, legend('Heunsistemi','ode45') subplot(2,2,2), hold on plot(t,Y(:,2),'r*'); legend('Heunsistemi','ode45') subplot(2,2,3), hold on plot(t,Y(:,3),'ro'); legend('Heunsistemi','ode45') subplot(2,2,4), hold on plot(t,Y(:,4),'r*'); legend('Heunsistemi','ode45') pause

% confronti con ode45 figure(3) subplot(2,2,1), hold on plot(t,Y(:,1),'ro');, legend('Rksistemi','ode45') subplot(2,2,2), hold on plot(t,Y(:,2),'r*'); legend('Rksistemi','ode45') subplot(2,2,3), hold on plot(t,Y(:,3),'ro');, legend('Rksistemi','ode45') subplot(2,2,4), hold on plot(t,Y(:,4),'r*'); legend('Rksistemi','ode45')

Dal command window >>mainSpacecraft

r cos()

r si

n(

)

r cos()

r si

n(

)

r cos()

r si

n(

)

r cos()

r si

n(

)

r cos()

r cos()

r cos()

r cos()

r si

n(

)

r si

n(

)

r si

n(

)

r si

n(

)

62

>> format long >> disp(T_RK(1:2,:)‘) 1.0e+06 * 0 7.150140000000000 0.000050000000000 7.148242144963171 0.000100000000000 7.142550735619172 0.000150000000000 7.133072266896018 0.000200000000000 7.119817656773666 0.000250000000000 7.102802387088111 0.000300000000000 7.082046702972614 0.000350000000000 7.057575873260109 0.000400000000000 7.029420514844428 0.000450000000000 6.997616984673679 0.000500000000000 6.962207843721299 0.000550000000000 6.923242397938913 0.000600000000000 6.880777321823856 0.000650000000000 6.834877370809193 0.000700000000000 6.785616189170178 0.000750000000000 6.733077220490102 0.000800000000000 6.677354727873715 0.000850000000000 6.618554930948934 0.000900000000000 6.556797266139562 0.000950000000000 6.492215775570221 0.001000000000000 6.424960628083602 0.001050000000000 6.355199772962839 0.001100000000000 6.283120722753624 0.001150000000000 6.208932455704089 0.001200000000000 6.132867420354973

RT =6378:14Km

1000s < t < 1050s

Esercizio

63

• Definire una griglia sul dominio D: >> [x,y]=meshgrid(-2:.1:2,-2:.1:2);

La funzione [x,y]=meshgrid(xvett,yvett) genera una griglia di punti le cui coordinate sono contenute nelle matrici x e y. I valori degli elementi di x e y sono dati dai vettori di punti equidistanti xvett e yvett

Grafici

64

x = -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 y = -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 -0.5000 0 0 0 0 0 0 0 0 0 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000

65

Grafici

• Definire f(x) e valutarla nei punti della griglia: >> f=@(x,y)[x.*exp(-x.^2-y.^2)];

>> z=f(x,y);

L’istruzione è equivalente a

>> z=x.*exp(-x.^2-y.^2);

66

Grafici

z =

-0.0007 -0.0029 -0.0067 -0.0071 0 0.0071 0.0067 0.0029 0.0007 -0.0039 -0.0167 -0.0388 -0.0410 0 0.0410 0.0388 0.0167 0.0039 -0.0135 -0.0582 -0.1353 -0.1433 0 0.1433 0.1353 0.0582 0.0135 -0.0285 -0.1231 -0.2865 -0.3033 0 0.3033 0.2865 0.1231 0.0285 -0.0366 -0.1581 -0.3679 -0.3894 0 0.3894 0.3679 0.1581 0.0366 -0.0285 -0.1231 -0.2865 -0.3033 0 0.3033 0.2865 0.1231 0.0285 -0.0135 -0.0582 -0.1353 -0.1433 0 0.1433 0.1353 0.0582 0.0135 -0.0039 -0.0167 -0.0388 -0.0410 0 0.0410 0.0388 0.0167 0.0039 -0.0007 -0.0029 -0.0067 -0.0071 0 0.0071 0.0067 0.0029 0.0007

67

Grafici

• Tracciare il grafico di z >> surf(x,y,z);

>> colorbar

Il comando surf(x,y,z) disegna la superficie z=f(x,y) valutata nei punti della griglia le cui coordinate sono date dagli elementi corrispondenti delle matrici x e y Colorbar mostra la barra dei colori che distinguono i valori assunti da z 68

Grafici

Il comando contour(x,y,z) disegna le linee di livello, cioè le curve dei punti in cui la superficie assume un valore fissato costante

69

Grafici

Il comando mesh(x,y,z) si può usare in alternativa al comando surf I comandi meshc(x,y,z) e surfc(x,y,z) disegnano contemporaneamente la superficie e le linee di livello

>> meshc(x,y,z)

70

Grafici

Il grafico si può completare aggiungendo etichette e titoli, esattamente come accade per il caso 1D Altre funzioni utili sono view cambia l’orientamento del grafico (punto di vista) colormap cambia il colore del grafico shading cambia l’ombreggiatura del grafico Usare lo help per stabilire quando usare le funzioni pcolor e plot3

71

Grafici

72

Grafici: esempi

73

Grafici: esempi

74

Grafici: esempi

75

Grafici: esempi

76

Esempio: f(x; y) = x2+y2-3, g(x,y) = xy-1

77

Esempio: f(x; y) = x2+y2-3, g(x,y) = xy-1

78

Esercizio

79

80

Sistema di equazioni non lineari

𝑓1 𝑥1, 𝑥2, 𝑥3 = 0

𝑓2 𝑥1, 𝑥2, 𝑥3 = 0

𝑓3 𝑥1, 𝑥2, 𝑥3 = 0

𝐶

1 + 𝑒 𝑠𝑖𝑛(𝛼 − 30)− 6870 = 0

𝐶

1 + 𝑒 𝑠𝑖𝑛(𝛼)− 6728 = 0

𝐶

1 + 𝑒 𝑠𝑖𝑛(𝛼 + 30)− 6615 = 0

𝐽𝐹 𝐶, 𝛼, 𝑒 =

1

1 + 𝑒𝑠𝑖𝑛(𝛼 − 30)−

𝐶 𝑒 cos (𝛼 − 30)

(1 + 𝑒𝑠𝑖𝑛(𝛼 − 30))2−

𝐶 sin (𝛼 − 30)

(1 + 𝑒𝑠𝑖𝑛(𝛼 − 30))2

1

1 + 𝑒𝑠𝑖𝑛(𝛼)−

𝐶 𝑒 cos (𝛼)

(1 + 𝑒𝑠𝑖𝑛(𝛼))2−

𝐶 sin (𝛼)

(1 + 𝑒𝑠𝑖𝑛(𝛼))2

1

1 + 𝑒𝑠𝑖𝑛(𝛼 + 30)−

𝐶 𝑒 cos (𝛼 + 30)

(1 + 𝑒𝑠𝑖𝑛(𝛼 + 30))2−

𝐶 sin (𝛼 + 30)

(1 + 𝑒𝑠𝑖𝑛(𝛼 + 30))2

Matrice jacobiana

F(X)=0

𝐶0, 𝛼0, 𝑒0 𝑎𝑝𝑝𝑟𝑜𝑠𝑠𝑖𝑚𝑎𝑧𝑖𝑜𝑛𝑒 𝑖𝑛𝑖𝑧𝑖𝑎𝑙𝑒

𝐽𝐹 𝐶𝑘 , 𝛼𝑘 , 𝑒𝑘 𝑉(𝑘) = 𝐹 𝐶𝑘 , 𝛼𝑘 , 𝑒𝑘 , 𝑘 = 0,1,2, …

𝑋(𝑘+1) = 𝑋(𝑘) − 𝑉 𝑘

Algoritmo

81

Per risolvere il problema con il metodo di Newton-Raphson è necessario: - Separare le radici in modo grafico

- Per ogni intervallo contenente una radice compatibile con il

problema fisico, verificare che il determinante della matrice JF non si annulla (Oss: det( JF) = 0 è un’equazione non lineare)

- Eseguire le iterazioni del metodo e determinare le incognite C, , e

Osservazione

82

%eliminando la dipendenza da C, riduciamo a due le equazioni del sistema

point_alfa=linspace(-pi/2,pi/2,50);

point_ecc = linspace(0,1,50);

[alfa,e]=meshgrid(point_alfa,point_ecc);

z1= 6870*(1+e.*sin(alfa-pi/6))./(1+e.*sin(alfa))-6728;

z2= 6870*(1+e.*sin(alfa-pi/6))./(1+e.*sin(alfa+pi/6))-6615;

figure,

subplot(2,2,1),meshc(alfa,e,z1), title('z1=z1(alfa,e)'), xlabel('alfa'), ylabel('e')

subplot(2,2,2),meshc(alfa,e,z2), title('z2=z2(alfa,e)'), xlabel('alfa'), ylabel('e')

subplot(2,2,3),contour(alfa,e,z1,[0 0],'g','linewidth',1),%linee di livello

z1(x,y)=0

hold on,contour(alfa,e,z2,[0 0],'r','linewidth',1) %linee di livello

z2(x,y)=0

xlabel('alfa'), ylabel('e'),legend('z1=0','z2=0'),title('Intersezioni')

Intervallo di separazione 0.03e0.05 /12/6

83

% script per la soluzione di un sistema non lineare di tre equazioni in tre % incognite (Problema orbita satellite script_minimoR.m) format long; % funzioni del sistema e derivate parziali rispetto alle incognite f1 = @(C,alfa,e)(C./(1+e.*sin(alfa-pi/6))-6870); f1x = @(C,alfa,e)(1./(1+e.*sin(alfa-pi/6))); f1y = @(C,alfa,e)(-C.*e.*cos(alfa-pi/6)./((1+e.*sin(alfa-pi/6)).^2)); f1z = @(C,alfa,e)(-C.*sin(alfa-pi/6)./((1+e.*sin(alfa-pi/6)).^2)); f2 = @(C,alfa,e)(C./(1+e.*sin(alfa))-6728); f2x = @(C,alfa,e)(1./(1+e.*sin(alfa))); f2y = @(C,alfa,e)(-C.*e.*cos(alfa)./((1+e.*sin(alfa)).^2)); f2z = @(C,alfa,e)(-C.*sin(alfa)./((1+e.*sin(alfa)).^2)); f3 = @(C,alfa,e)(C./(1+e.*sin(alfa+pi/6))-6615); f3x = @(C,alfa,e)(1./(1+e.*sin(alfa+pi/6))); f3y = @(C,alfa,e)(-C.*e.*cos(alfa+pi/6)./((1+e.*sin(alfa+pi/6)).^2)); f3z = @(C,alfa,e)(-C.*sin(alfa+pi/6)./((1+e.*sin(alfa+pi/6)).^2));

84

% punto iniziale xk = 6870; yk = pi/6; zk=0.03; % inizializzazione variabili iter = 0; errC = 10.; erralfa = 10.; erre = 10.; errf1 = abs(f1(xk,yk,zk)); errf2 = abs(f2(xk,yk,zk)); errf3=abs(f3(xk,yk,zk)); [iter xk errC errf1 yk erralfa errf2 zk erre errf3] % iterazioni del metodo di Newton (numero di iterazioni prefefinite) for iter= 1:10 xv = xk; yv = yk; zv = zk; JFk = [f1x(xv,yv,zv),f1y(xv,yv,zv),f1z(xv,yv,zv);f2x(xv,yv,zv),f2y(xv,yv,zv),f2z(xv,yv,zv);… f3x(xv,yv,zv),f3y(xv,yv,zv),f3z(xv,yv,zv)]; % matrice jacobiana [iter det(JFk)] Bk = [f1(xv,yv,zv);f2(xv,yv,zv);f3(xv,yv,zv)]; Vk = JFk\Bk xk = xv-Vk(1);, yk = yv-Vk(2);, zk = zv-Vk(3); errC = abs(xk-xv); erralfa = abs(yk-yv); erre = abs(zk-zv); errf1 = abs(f1(xk,yk,zk)); errf2 = abs(f2(xk,yk,zk)); errf3 =abs(f3(xk,yk,zk)); [iter xk errC errf1 yk erralfa errf2 zk erre errf3], pause end disp(‘soluzione’) [xk yk zk]

Modificare usando il comando input

85

Dal command window >> script_minimoR >> format long >> [xk yk zk] ans = 1.0e+03 * 6.819293793207854 0.000340783997950 0.000040598959058 >> R=xk/(1+zk) R = 6.553239107006842e+03 >> theta=pi/2-yk theta = 1.230012328844540

Esercizio Scrivere la funzione Matlab SOR_opt.m che implementi il metodo iterativo S.O.R.. La funzione deve ricevere in input la matrice A del sistema, il vettore b dei termini noti, il vettore dell’approssimazione iniziale X0, l’accuratezza eps richiesta alla approssimazione della soluzione prodotta e il numero massimo di iterazioni consentite max_iter.

Se il parametro di rilassamento non è dato in input, la funzione deve assegnare alla variabile omega il valore del parametro ottimo.

La funzione deve restituire in output, la approssimazione X della soluzione del sistema, il vettore ERR dell’errore massimo commesso ad ogni iterazione e il numero di iterazioni eseguite iter.

Se richiesto, la funzione deve restituire come variabile output il parametro di rilassamento omega usato.

Esercizio Scrivere la funzione Matlab SOR_opt.m che implementi il metodo iterativo S.O.R.. La funzione deve ricevere in input la matrice A del sistema, il vettore b dei termini noti, il vettore dell’approssimazione iniziale X0, l’accuratezza eps richiesta alla approssimazione della soluzione prodotta e il numero massimo di iterazioni consentite max_iter.

Se il parametro di rilassamento non è dato in input, la funzione deve assegnare alla variabile omega il valore del parametro

ottimo. La funzione deve rifiutare valori di omega (0,2)

La funzione deve restituire in output, la approssimazione X della soluzione del sistema, il vettore ERR dell’errore massimo commesso ad ogni iterazione e il numero di iterazioni eseguite iter.

Se richiesto, la funzione deve restituire come variabile output il parametro di rilassamento omega usato.

???

Strutture dati Cell permette di collezionare in un’unica variabile: •oggetti di vario tipo (vettori, matrici, variabili numeriche o logiche, caratteri o stringhe) •variabili dello stesso tipo ma di dimensione diversa

Variabili di tipo cell si definiscono tra parentesi graffe. Le componenti si elencano una dopo l’altra e separate da virgole C = {componente1, componente2, …. , componenteN}; C{i} estrae la i-esima componente di C C{i}(k) estrae (se esiste) il k-simo elemento della i-esima componente di C

Strutture dati Esempio: >> C = ,‘ciao’, *4 3 1 2+,3,*1 7 2; 0 5 8; 1 0 9+-;

>> whos Name Size Bytes Class C 1x4 488 cell array >>C{2} % estrae il secondo elemento della variabile di tipo cell C ans = 4 3 1 2 >>C{4}(2,3) %estrae l’elemento con indice di riga 2 e indice di colonna 3 della quarta componente di C ans = 8

definisce una variabile cell composta da 4 elementi: • una stringa • un vettore • un numero • una matrice

Strutture dati

Nota: una variabile di tipo Cell contiene copie di variabili e non puntatori a variabili. Quindi se C ={A,B}, con A, B due generiche variabili, il contenuto di C non cambia se in seguito A e B sono modificate Esempio: >> A = [2 5 3 6]; >> B = [1 5; 8 9;3 6]; >> C = {A,B}; >> B = 1; >> C{2} ans = 1 5 8 9 3 6

Strutture dati C può essere inizializzata con il comando C = cell(m,n) : C è composta da mxn matrici vuote Celldisp(C): visualizza il contenuto di una cell Cellplot(C): disegna il contenuto di una cell Le celle possono essere annidate, cioè un elemento di una cell può essere esso stesso una cell. Esempio: C = cell(1,4); A = [1 2 3; 4 5 6];, v = [1 4 6]; C{1} = {A,v}; C{2} = 4; C,4- = ‘stringa’; Per estrarre il contenuto della matrice A si digita il comando C{1}{1} gli indici sono ordinati da sinistra verso destra, dalla cell più esterna a quella più interna

Alcune funzioni

varargin = variabile di tipo cell contenente le variabili di input di una funzione. Si usa quando gli input di una funzione possono variare (esistono parametri di default) Oss: l’uso di varargin richiede un ordine preciso delle variabili di input. L’assegnazione dei valori di default deve essere compresa nel primo blocco di istruzioni della funzione combinata con il comando nargin. nargin = numero delle variabili di input

Alcune funzioni Esempio: function [R] = confronta(x,varargin) % la funzione stabilisce se un vettore x ha elementi maggiori di % un valore fissato T. Se T non viene dato in input, viene assegnato il % valore di default 256. Lo output è la variabile logica R che vale 1 se % esiste almeno un elemento in x che risulta maggiore di T, 0 altrimenti. if nargin == 0 error(‘Attenzione: e’’ necessario introdurre un vettore’) elseif nargin == 1 T = 256; elseif nargin == 2 T = varargin{1}; else error(‘Troppe variabili di input!!!’) end ind = find(x>T); if isempty(ind), R = 0, else R=1; end

Alcune funzioni OSS: per lo output si usa varargout combinata con nargout (numero delle variabili di output) Esempio: function [R,varargout] = confronta(x,varargin) % la funzione stabilisce se un vettore x ha elementi maggiori di % un valore fissato T. Se T non viene dato in input, viene assegnato il valore di % default 256. Lo output è la variabile logica R che vale 1 se esiste almeno % un elemento > T, 0 altrimenti. La seconda variabile di output, se richiesta, % è il vettore contenente le posizioni degli elementi di % x > T ……… ……… ind = find(x>T); if isempty(ind), R = 0, else R=1; end if nargout == 2 varargout{1} = ind; elseif nargout>2 error(‘Troppe variabili di output’) end

Come funzione precedente

95

Esercizio 1

• Scrivere la funzione Matlab num_cond.m che legga in input una matrice quadrata A e la variabile numerica tipo_norma che può essere uguale a 1, 2 o inf. La funzione restituisca in output il numero di condizionamento rispetto alla norma indicata nella variabile tipo_norma. Se quest’ ultima non viene data in input, la funzione deve calcolare il numero di condizionamento rispetto alla norma infinito.

Si usi la funzione num_cond.m per calcolare il numero di condizionamento delle seguenti matrici

199.0

01.11 A

199

1011 A

40

81

32

A

96

Soluzione

• Determinare la lista delle variabili di input e di output

Variabili di input: - matrice A di cui calcolare in numero di condizionamento - tipo di norma (tipo_norma) rispetto al quale calcolare il condizionamento. tipo_norma può non essere data in input

Variabili di output: - numero di condizionamento K Controllo sugli input: - controllo sul numero degli input - la matrice deve essere quadrata - il tipo di norma deve essere un numero pari a 1, 2, o inf

97

• Aprire la finestra di Editor

• Salvare il file come num_cond.m nella directory di lavoro • Scrivere la seguente istruzione che contiene la lista degli output e degli input della funzione function [K] = num_cond(A,varargin) • Scrivere lo help della funzione % function [K] = num_cond(A,varargin) % calcola il numero di condizionamento della matrice quadrata % A rispetto alla norma indicata in tipo_norma. % Se tipo_norma non è data in input, il condizionamento è calcolato % rispetto alla norma infinito. tipo_norma non può assumere valori % diversi da ‘1’, ‘2’ o ’inf’

La funzione prevede parametri di input opzionali!

98

• Controllo del numero delle variabili di input e controllo degli input

if nargin == 0 % non sono stati assegnati input error('E'' necessario definire una matrice!!!') elseif nargin == 1 % è stato passato un solo input tipo_norma = inf; elseif nargin == 2

tipo_norma = varargin{1}; if (tipo_norma ~= 1) & (tipo_norma ~= 2) & (tipo_norma ~= inf) error('la norma deve essere 1 2 o inf') end

else error('troppe variabili di input!!!')

end

Nota: poiché è stato assegnato un solo

input, allora necessariamente è la matrice (o comunque

deve essere assegnato alla variabile che corrisponde alla

matrice). In questo caso è necessario

assegnare il valore di default alla variabile

tipo_norma, necessaria per il calcolo del

numero di condizionamento

Nota: poiché sono stati assegnati due input,

allora necessariamente il primo è la matrice

mentre il secondo deve essere assegnato alla variabile tipo_norma,

99

• Controllo del numero delle variabili di input

if nargin == 0 % non sono stati assegnati input error('E'' necessario definire una matrice!!!') elseif nargin == 1 % è stato passato un solo input tipo_norma = inf; elseif nargin == 2

tipo_norma = varargin{1}; if (tipo_norma ~= 1) & (tipo_norma ~= 2) & (tipo_norma ~= inf) error('la norma deve essere 1 2 o inf') end

else error('troppe variabili di input!!!') end

Nota: Lo input da assegnare a tipo_norma

è il primo elemento della cell varargin, cioè quella contenente tutte

le variabili di input opzionali con cui è stata chiamata la

funzione

Nota: poiché sono stati assegnati più di due

input, si da il messaggio di errore

100

• Scrivere la lista di istruzioni per il controllo degli input dim = size(A); if dim(1)~=dim(2) error('La matrice deve essere quadrata') end if det(A)==0 error('La matrice e’’ singolare') end

• Scrivere la lista di istruzioni per il calcolo del numero di condizionamento rispetto al tipo di norma scelto

IA = inv(A); % matrice inversa di A NA = norm(A,tipo_norma); % norma di A NIA = norm(IA,tipo_norma); % norma della matrice inversa K = NA*NIA; % numero di condizionamento

• Salvare il file

Nota: il controllo su tipo_noma è stato fatto nel caso in cui nargin=2 quindi non va ripetuto

101

Dalla finestra del Command Window >> A = [1 1.01;0.99 1]; >> KA_1 = num_cond(A,1); >> disp(KA_1) 4.0401e+004 >> KA = num_cond(A); >> disp(KA) 4.0401e+004 >> KA = num_cond(A,-1); ??? Error using ==> num_cond at 18 la norma deve essere 1 2 o inf >> A = [2 3; -1 8; 0 4]; >>KA_inf =num_cond(A,inf); ??? Error using ==> num_cond at 10 La matrice deve essere quadrata

Non è stato specificato il tipo di norma

il valore di input

assegnato per tipo di norma

non è accettabile

102

Esercizio 2

• Modificare la funzione Matlab num_cond.m e chiamarla num_cond2.m. La

funzione deve prevedere come ulteriori variabili di output, se richieste, la matrice inversa di A e la norma di A calcolata rispetto a quanto indicato in tipo_norma.

103

• Aprire la finestra di Editor

• Salvare il file come num_cond2.m nella directory di lavoro • Scrivere la seguente istruzione che contiene la lista degli output e degli input della funzione function [K,varargout] = num_cond2(A,varargin) • Scrivere lo help della funzione % function [K,varargout] = num_cond2(A,varargin) % calcola il numero di condizionamento K della matrice quadrata % A rispetto alla norma indicata in tipo_norma. % Se tipo_norma non è data in input, il condizionamento è calcolato % rispetto alla norma infinito. tipo_norma non può assumere valori % diversi da ‘1’, ‘2’ o ’inf’ % Se si chiedono due output, allora il primo è K e il secondo è la % matrice inversa di A % Se si chiedono tre output, il primo è K, il secondo è la matrice % inversa di A, il terzo è la norma di A

La funzione prevede parametri opzionali sia di input che di

output

104

if nargin == 0 % non sono stati assegnati input error('E'' necessario definire una matrice!!!') elseif nargin == 1 % è stato passato un solo input tipo_norma = inf; elseif nargin == 2 % sono stati specificati due input

tipo_norma = varargin{1}; if (tipo_norma ~= 1) & (tipo_norma ~= 2) & (tipo_norma ~= inf) error('la norma deve essere 1 2 o inf') end

else % sono stati assegnati input più di due input error('troppe variabili di input!!!')

end

• Lista di istruzioni per il controllo degli input dim = size(A); if dim(1)~=dim(2) error('La matrice deve essere quadrata') end if det(A)==0 error('La matrice e’’ singolare') end

• Controllo del numero delle variabili di input

105

• Lista di istruzioni per il calcolo del numero di condizionamento rispetto al tipo di norma scelto

IA = inv(A); % matrice inversa di A NA = norm(A,tipo_norma); % norma di A NIA = norm(IA,tipo_norma); % norma della matrice inversa K = NA*NIA; % numero di condizionamento

• Controllo e assegnazioni output

if nargout == 2 % sono stati richiesti due output varargout{1} = IA; elseif nargout == 3 %sono stati richiesti tre output varargout{1} = IA; varargout{2} = NA; elseif nargout >3% sono stati richiesti troppi output

error('troppe variabili di output!!!') end

Nota: poiché sono stati richiesti due

output, allora necessariamente il primo è K mentre il secondo è l’inversa di A, che quindi diventa il primo elemento della

variabile di tipo cell varargout, che

contiene tutti gli output opzionali

richiesti

106

• Lista di istruzioni per il calcolo del numero di condizionamento rispetto al tipo di norma scelto

IA = inv(A); % matrice inversa di A NA = norm(A,tipo_norma); % norma di A NIA = norm(IA,tipo_norma); % norma della matrice inversa K = NA*NIA; % numero di condizionamento

• Controllo e assegnazioni output

if nargout == 2 % sono stati richiesti due output varargout{1} = IA; elseif nargout == 3 %sono stati richiesti tre output varargout{1} = IA; varargout{2} = NA; elseif nargout >3% sono stati richiesti troppi output

error('troppe variabili di output!!!') end

Nota: poiché sono stati richiesti tre

output, allora necessariamente il

primo è K, il secondo è l’inversa di A, che diventa il primo elemento della variabile di tipo cell varargout, mentre il terzo è la norma di A che diventa il

secondo elemento di varargout

107

Dalla finestra del Command Window >> A = [1 1.01;0.99 1]; >> KA_1 = num_cond2(A,1); >> disp(KA_1) 4.0401e+004 >> [KA,B] = num_cond2(A); >> disp(KA) 4.0401e+004 >> disp(B) 1.0e+004 * 1.0000 -1.0100 -0.9900 1.0000 >> [KA,B,nA] = num_cond2(A); >> disp(nA) 2.0100

Non è stato specificato il tipo di norma

E’ richiesto un terzo output

E’ richiesto un secondo output

Esercizio Tornando all'esercizio proposto relativo al metodo SOR function [X,ERR,iter,varargout] = SOR_opt(A,b,X0,eps,max_iter,varargin) % function [X, ERR, iter] = SOR_opt(A,b,X0,eps,max_iter,varargin) % Risolve un sistema lineare con metodo di sovra-rilassamento SOR % eventualmente usando il parametro ottimo qualora questo non sia % un parametro di input % INPUT: % A = matrice dei coefficienti del sistema % B = vettore dei termini noti % X0 = vettore dell'approssimazione iniziale % eps = accuratezza della soluzione % max_iter = numero massimo di iterazioni consentite % varargin{1} = omega, parametro di rilassamento del metodo % OUTPUT: % X = vettore soluzione % ERR = vettore dell'errore massimo per ogni iterazione % iter = numero di iterazioni eseguite % varargout{1} = valore del parametro di rilassamento

% controllo e assegnazione valori di default alle variabili di input if nargin == 5 L = tril(A); Minv = inv(L); U = triu(A,1); C_GS = -Minv*(U); rho_GS = max(abs(eig(C_GS))); if (rho_GS >= 1) error('Attenzione: ==>> rho di Gauss Seidel > 1') else omega = 2/(1+sqrt(1-rho_GS)); end elseif nargin == 6 omega = varargin{1}; if (omega<=0 | omega >= 2) error(‘il parametro di rilassamento omega non e'' in (0,2)') end df = def_pos(A);

Modificare includendo un controllo per nargin<5 e sul determinante della matrice

if df == 0 % Costruzione matrice di iterazione L = tril(A,-1); D = diag(diag(A)); Minv = inv(L+D/omega); U = triu(A,1); C_sor = -Minv*(U-(1-omega)*D/omega); % Verifica C.N.S. di convergenza per i metodi iterativi rhoCsor = max(abs(eig(C_sor))) if (rhoCsor >= 1) error('Attenzione: ==>> rho > 1, il metodo non converge') end end elseif nargin > 6 error('troppe variabili di input') end

% controllo e assegnazione valori di default alle variabili di output if nargout == 4 varargout{1} = omega; elseif nargout > 4 error('troppe variabili di output') end % Calcolo delle dimensioni della matrice dimA = size(A);, n = dimA(1);, X0 = X0(:)'; % Ciclo iterativo err = 100;, iter = 0;, ERR = [];, X = X0; tic, while (err>eps & iter<= max_iter) for i = 1:n V(i)=(-sum(A(i,1:i-1).*X(1:i-1))-sum(A(i,i+1:n).*X0(i+1:n))+b(i))/A(i,i); X(i)=omega * V(i) + (1 - omega)*X0(i); end err = norm(X-X0,inf); ERR = [ERR err];, X0 = X;, iter = iter + 1; end, toc, if (iter > max_iter) & (err>eps) fprintf('accuratezza richiesta non raggiunta dopo %7d iterazioni',max_iter) end

function [df] = def_pos(A) %function [df] = def_pos(A) % stabilisce se la matrice quadrata A è simmetrica e definita positiva usando il % criterio di Sylvester % INPUT % A = matrice quadrata % OUTPUT % df = variabile logica. df = 1 se A è simmetrica e definita positiva; 0 altrimenti [m,n] = size(A); % controlla se A è una matrice quadrata if m ~=n, error('la matrice A deve essere quadrata!!!'), end S = A==A‘; if all(S) i = 1; while (i<=m) & (det(A(1:i,1:i))>0) % controllo prima la dimensione e poi il determinante!!! i = i+1; end if i==m+1 df = 1;, disp('la matrice e'' simmetrica e definita positiva') else df = 0;, disp('la matrice e'' simmetrica ma non e'' definita positiva') end else disp('la matrice non e'' simmetrica'), df = 0; end

Switch

Switch-case-Otherwise:

switch nome_variabile (o espressione)

case valore1

1° blocco di istruzioni

case valore2

2° blocco di istruzioni

……….

otherwise

ultimo blocco di istruzioni

end

Se il valore di nome_variabile è valore1 viene eseguito il 1° blocco di

istruzioni;

Se è valore2, il 2° blocco di istruzioni, e così via,

altrimenti esegue l’ultimo blocco di istruzioni

Switch

Esempio: a = rem(b,3); % resto della divisione per 3

switch a

case 0

c = b/3;

case 1

c = (b-a)/3;

otherwise

c = b-a/2;

end

118

Esercizio

• Modificare la funzione Matlab num_cond2.m e chiamarla num_cond3.m

usando opportunamente l’istruzione switch/case/otherwise/end

119

function [K,varargout] = num_cond3(A,varargin) % function [K,varargout] = num_cond3(A,varargin) % calcola il numero di condizionamento K della matrice quadrata % A rispetto alla norma indicata in tipo_norma. % Se tipo_norma non è data in input, il condizionamento è calcolato % rispetto alla norma infinito. tipo_norma non può assumere valori % diversi da ‘1’, ‘2’ o ’inf’ % Se si chiedono due output, allora il primo è K e il secondo è la % matrice inversa di A % Se si chiedono tre output, il primo è K, il secondo è la matrice % inversa di A, il terzo è la norma di A

A parte il nome, in questa prima parte non

cambia nulla rispetto alla funzione num_cond2

120

switch nargin case 0 % non sono assegnati input error('E'' necessario definire una matrice!!!') case 1 % è assegnato un solo input, cioè la matrice tipo_norma = inf; case 2 % sono assegnati due input, cioè la matrice e tipo_norma tipo_norma = varargin{1};

if (tipo_norma ~= 1) & (tipo_norma ~= 2) & (tipo_norma ~= inf) error('la norma deve essere 1 2 o inf') end

otherwise error('troppe variabili di input!!!')

end

• Controllo del numero delle variabili di input

L’istruzione condizionale annidata

if elseif elseif else end usata in num_cond2 è

sostituita con l’istruzione switch sulla variabile

nargin

121

• Lista di istruzioni per il controllo degli input

dim = size(A); if dim(1)~=dim(2) error('La matrice deve essere quadrata') end if det(A)==0 error('La matrice e’’ singolare') end

• Lista di istruzioni per il calcolo del numero di condizionamento rispetto al tipo di norma scelto

IA = inv(A); % matrice inversa di A NA = norm(A,tipo_norma); % norma di A NIA = norm(IA,tipo_norma); % norma della matrice inversa K = NA*NIA; % numero di condizionamento

non cambia nulla rispetto alla funzione num_cond2

122

• Controllo e assegnazioni output

switch nargout case 2 % sono stati richiesti due output varargout{1} = IA; case 3 %sono stati richiesti tre output varargout{1} = IA; varargout{2} = NA; otherwise % il numero delle variabili di output richieste è troppo % alto o troppo basso

error('controllare il numero delle variabili di output!!!') end

• L’istruzione condizionale annidata if elseif else end usata in num_cond2 è sostituita con l’istruzione switch sulla variabile nargout

• otherwise corrisponde sia ad un numero di output troppo grande rispetto a quello previsto, sia ad un numero troppo piccolo!!!

123

Dalla finestra del Command Window >> A = [1 1.01;0.99 1]; >> KA_1 = num_cond3(A,1); >> disp(KA_1) 4.0401e+004 >> [KA,B] = num_cond3(A); >> disp(KA) 4.0401e+004 >> disp(B) 1.0e+004 * 1.0000 -1.0100 -0.9900 1.0000 >> [KA,B,nA] = num_cond3(A); >> disp(nA) 2.0100

non cambia nulla rispetto alla funzione num_cond2

Esercizio Modificare la funzione SOR_opt.m assegnando valori di default alle variabili X0, eps, max_iter e omega nel caso in cui queste non siano date in input. La funzione restituisca in output la soluzione e, solo se richiesto, l'errore ad ogni iterazione, il numero di iterazioni eseguite e il valore del parametro di rilassamento. Si usi opportunamente il comando switch-case-otherwise sia sulla variabile nargin che su naragout.

function [X,varargout] = SOR_opt_varargin(A,b,varargin) ……………………………………………..…………………………………………….. ……………………………………………..…………………………………………….. if nargin<2 error('I dati di input sono incompleti') end switch nargin case 2 X0 = zeros(length(b),1)’; eps = 0.5*10^-5; max_iter = 100; L = tril(A);, Minv = inv(L);, U = triu(A,1);, C_GS = -Minv*(U); rho_GS = max(abs(eig(C_GS))); if (rho_GS >= 1) error('Attenzione: ==>> rho di Gauss Seidel > 1') else omega = 2/(1+sqrt(1-rho_GS)); end

case 3 X0 = varargin{1}; eps = 0.5*10^-5; max_iter = 100; L = tril(A);, Minv = inv(L);, U = triu(A,1);, C_GS = -Minv*(U); rho_GS = max(abs(eig(C_GS))); if (rho_GS >= 1) error('Attenzione: ==>> rho di Gauss Seidel > 1') else omega = 2/(1+sqrt(1-rho_GS)); end case 4 X0 = varargin{1}; eps = varargin{2}; max_iter = 100; L = tril(A);, Minv = inv(L);, U = triu(A,1);, C_GS = -Minv*(U); rho_GS = max(abs(eig(C_GS))); if (rho_GS >= 1), error('Attenzione: ==>> rho di Gauss Seidel > 1') else, omega = 2/(1+sqrt(1-rho_GS));, end

case 5 X0 = varargin{1}; eps = varargin{2}; max_iter = varargin{3}; L = tril(A);, Minv = inv(L);, U = triu(A,1);, C_GS = -Minv*(U); rho_GS = max(abs(eig(C_GS))); if (rho_GS >= 1) error('Attenzione: ==>> rho di Gauss Seidel > 1') else omega = 2/(1+sqrt(1-rho_GS)); end case 6 X0 = varargin{1}; eps = varargin{2}; max_iter = varargin{3}; omega = varargin{4}; if (omega<=0 | omega >= 2) error('il metodo SOR non converge: omega non e'' in (0,2)') end

if df == 0 % Costruzione matrice di iterazione L = tril(A,-1); D = diag(diag(A)); Minv = inv(L+D/omega); U = triu(A,1); C_sor = -Minv*(U-(1-omega)*D/omega); % Verifica C.N.S. di convergenza per i metodi iterativi rhoCsor = max(abs(eig(C_sor))) if (rhoCsor >= 1) error('Attenzione: ==>> rho > 1, il metodo non converge') end end otherwise error('Troppe variabili di input!!!') end … …

Implementazione del metodo

if nargout >1 switch nargout case 2 varargout{1} = ERR; case 3 varargout{1} = ERR; varargout{2} = iter; case 4 varargout{1} = ERR; varargout{2} = iter; varargout{3} = omega; otherwise error('Troppe variabili di input!!!') end end