Metodi Numerici con elementi di Programmazione

Post on 30-Jul-2022

5 views 0 download

Transcript of Metodi Numerici con elementi di Programmazione

Metodi Numerici con elementi di Programmazione

A.A. 2019-2020

Introduzione al MatLab

(4)

1

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

Durata di un programma Per confrontare due programmi che risolvono lo stesso problema è utile misurare il tempo di CPU impiegato per eseguirli. La funzione cputime misura il tempo di CPU usato da quando è stata aperta la sessione di Matlab Per misurare il tempo richiesto dall’ esecuzione di un programma, quindi, basta valutare la differenza tra il cputime valutato prima e dopo l’esecuzione del programma In alternativa si possono usare in modo combinato i comandi tic e toc tic si scrive all’ inizio del programma, o immediatamente prima della prima istruzione da cui si vuole iniziare a valutare il tempo di esecuzione toc si scrive alla fine del programma, o immediatamente dopo l’ultima istruzione della parte di programma di cui si vuole valutare il tempo di esecuzione. toc contiene il tempo di esecuzione In questa caso il tempo di esecuzione viene stampato sulla finestra dei comandi

Esercizio Scrivere lo script Matlab tempo_prodotto.m che crei una matrice quadrata A di dimensione n e un vettore v di lunghezza n e ne calcoli il prodotto matriciale. Gli elementi della matrice e del vettore sono numeri casuali; la dimensione n deve essere fatta variare da n=150 a n=9000 usando un passo di 150. Misurare e graficare opportunamente il tempo di esecuzione del prodotto A*v

Esercizio % Lo script tempo_prodotto.m valuta il tempo di calcolo richiesto dalla % valutazione del prodotto matrice per vettore al variare della dimensione % della matrice n_in= 150; , n_fin = 9000; p = 150; A = rand(n); V = rand(n,1); T = []; for k=n_in:p:n_fin B = A(1:k,1:k); w = V(1:k); ti = cputime; P = B*w; tf = cputime-ti; T = [T tf]; end figure, plot(n_in:p:n_fin,T), xlabel(‘dimensione della matrice’), ylabel(‘tempo’)

Esercizio % Lo script tempo_prodotto.m valuta il tempo di calcolo richiesto dalla % valutazione del prodotto matrice per vettore al variare della dimensione % della matrice n_in= 150; , n_fin = 9000; p = 150; A = rand(n); V = rand(n,1); T = []; for k=n_in:p:n_fin B = A(1:k,1:k); w = V(1:k); tic; P = B*w; toc; T = [T toc]; end figure, plot(n_in:p:n_fin,T), xlabel(‘dimensione della matrice’), ylabel(‘tempo’)

8

Scrivere lo script Matlab MEG_nopivot.m che implementa il metodo di eliminazione di Gauss senza pivoting e confrontare il risultato con quello prodotto usando l’operatore \

9

10

La matrice è regolare (det(A) = 1.499999999999999e+06) Inoltre, det(Ak) 0, k=1,…,6; infatti det(A1) = 11 det(A2) = 351 det(A3)=1962 det(A4)= 2520 det(A5)=44999.99999999999e+04 det(A6)=det(A) = 1.499999999999999e+06

Il vettore dei termini noti è b = [500 0 0 0 0 0]’

11

Si eseguono n-1 passi

% meg_nopivot % Risolve il sistema lineare Ax=b usando il metodo di % eliminazione di Gauss senza pivoting

clear all % cancella tutte le variabili in Workspace %lettura della matrice dei coefficienti e del vettore dei termini noti A = input(‘inserisci la matrice dei coefficienti del sistema A= ’); b = input(‘inserisci il vettore dei termini noti b = ‘’); [m,n] = size(A); % controllo sulla dimensione di A e b if (m˜=n) || ((m==n) && (n˜=length(b))) error(‘Attenzione alla dimensione di A e b!!!’); end % controllo dell’esistenza di un’unica soluzione if det(A) == 0 error(‘la matrice dei coefficienti e’’ singolare!!!’) end

Soluzione (script meg_nopivot.m)

12

13

%controllo dei minori principali di testa for i=1:n-1 if det(A(1:i,1:i)) == 0 error('uno dei minori principali di testa e'' singolare!!!') end end for k =1:n-1 % gestisce il numero di passi del metodo di Gauss for i =k+1:n m = A(i,k)/A(k,k); % calcolo dei moltiplicatori for j=k:n % modifica degli elementi sulla i-sima riga A(i,j) = A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); % modifica della i-sima riga componente di b end end

14

%controllo dei minori principali di testa for i=1:n-1 if det(A(1:i,1:i)) == 0 error('uno dei minori principali di testa e'' singolare!!!') end end for k =1:n-1 % gestisce il numero di passi del metodo di Gauss for i =k+1:n m = A(i,k)/A(k,k); % calcolo dei moltiplicatori for j=k:n % modifica degli elementi sulla i-sima riga A(i,j) = A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); % modifica della i-sima riga componente di b end end

Non è necessario

controllare che

A(k,k) sia diverso

da zero !!!

% metodo di sostituzione all’indietro % All’uscita dei cicli annidati precedenti la matrice A è triangolare superiore x(n) = b(n)/A(n,n); for i =n-1:-1:1 x(i) = (b(i)-sum(A(i,i+1:n).*x(i+1:n)))/A(i,i); end x = x’; % la soluzione deve essere un vettore colonna

fprintf('la soluzione del sistema e'' x= \n ')

fprintf('%+6.5f \n',x(:))

15

16

Dal Command Window >> MEG_nopivot inserisci la matrice del sistema A= [11 -5 0 0 0 -1; -20 41 -15 0 -6 0; 0 -3 7 -4 0 0; 0 0 -1 2 -1 0; 0 -3 0 -10 28 -15; -2 0 0 0 -15 47]; inserisci il vettore dei termini noti b = *500 0 0 0 0 0+’; la soluzione del sistema e' x= +70.00000 +52.00000 +40.00000 +31.00000 +22.00000 +10.00000 >> MEG_nopivot inserisci la matrice del sistema A= [11 -5 0 0 0 -1; 0 -3 0 -10 28 -15; -2 0 0 0 -15 47; 0 0 -1 2 -1 0; -20 41 -15 0 -6 0; 0 -3 7 -4 0 0] inserisci il vettore dei termini noti b = [500 0 0 0 0 0]' Error using MEG_nopivot (line 25) uno dei minori principali di testa e' singolare!!!

la quinta e la quarta equazione sono state scambiate con la seconda e la terza

17

Dal Command Window OSS: >> X=A\b X = 70.0000 52.0000 40.0000 31.0000 22.0000 10.0000 >> MEG_nopivot inserisci la matrice del sistema A= [11 -5 0 0 0 -1; -20 41 -15 0 -6 0; 0 -3 7 -4 0 0; 0 0 -1 2 -1 0; 0 -3 0 -10 28 -15; -2 0 0 0 -15 47]; inserisci il vettore dei termini noti b = [500 0 0 0 0 +’; ??? Error using ==> MEG_nopivot at 14 Attenzione alla dimensione di A e b!!!

%controllo dei minori principali di testa for i=1:n-1 if det(A(1:i,1:i)) == 0 disp(‘Attenzione: uno dei minori principali di testa e'' singolare!!!') break end end for k =1:n-1 % gestisce il numero di passi del metodo di Gauss for i =k+1:n if A(k,k) ˜=0 m = A(i,k)/A(k,k); % calcolo dei moltiplicatori for j=k:n % modifica degli elementi sulla i-sima riga A(i,j) = A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); % modifica della i-sima riga componente di b else error(‘E’’ necessario il pivoting’) end end

Osservazione 1: modifiche allo script

18

Non si interrompe l’esecuzione dello script ma solo del ciclo for. Per questo motivo è necessario controllare che A(k,k) sia diverso da zero quando si calcolano i moltiplicatori !!!

metodo di sostituzione all’indietro x(n) = b(n)/A(n,n); for i =n-1:-1:1 x(i) = (b(i)-sum(A(i,i+1:n).*x(i+1:n)))/A(i,i); (*) end Nota: 1. l’istruzione precedente è corretta perché x è generato come un

vettore riga 2. Il ciclo precedente è equivalente a for i =n-1:-1:1 somma = 0; for j=i+1:n somma = somma + A(i,j)*x(j); end x(i)=(b(i)-somma)/A(i,i); end 3. Stabilire se la seguente istruzione è equivalente a (*)

x(i) = (b(i)-A(i,i+1:n)*x(i+1:n)’)/A(i,i);

Osservazione 2: modifiche allo script

19

20

>> MEG_nopivot inserisci la matrice del sistema A= [-sqrt(3)/2 0 .5 0 0 0; -.5 0 -sqrt(3)/2 0 0 0; sqrt(3)/2 1 0 1 0 0; 0.5 0 0 0 1 0; 0 -1 -.5 0 0 0; 0 0 sqrt(3)/2 0 0 1] inserisci il vettore dei termini noti b = [0 500 0 0 0 0]' Error using MEG_nopivot (line 25) uno dei minori principali di testa e' singolare!!!

Dal Command Window

OSS: >> X=A\b X = -250.0000 216.5064 -433.0127 0 125.0000 375.0000

21

>> MEG_pivot inserisci la matrice del sistema A= [-sqrt(3)/2 0 .5 0 0 0; -.5 0 -sqrt(3)/2 0 0 0; sqrt(3)/2 1 0 1 0 0; 0.5 0 0 0 1 0; 0 -1 -.5 0 0 0; 0 0 sqrt(3)/2 0 0 1] inserisci il vettore dei termini noti b = [0 500 0 0 0 0]' la soluzione del sistema e' x= -250.00000 +216.50635 -433.01270 +0.00000 +125.00000 +375.00000

Dal Command Window

OSS: >> X=A\b X = -250.0000 216.5064 -433.0127 0 125.0000 375.0000

Alcune funzioni D=diag(v) costruisce una matrice diagonale con il vettore v sulla diagonale principale v=diag(D) estrae la diagonale principale della matrice D w=diag(D,k) se k > 0 estrae la k-esima diagonale superiore, se k < 0 estrae la k-esima diagonale inferiore T = triu(A) estrae la parte triangolare superiore di A. T è una matrice triangolare superiore T = triu(A,k) estrae gli elementi che sono al di sopra della k-esima diagonale di A. T = tril(A) estrae gli elementi che sono al di sotto della k-esima diagonale di A.

Alcune funzioni

R = rank(A) calcola il numero di righe linearmente indipendenti di A N = norm(A,n) calcola la norma n di A. Se n = 1, calcola la norma 1 n = 2 , calcola la norma 2 n = inf, calcola la norma infinito B=inv(A) calcola la matrice inversa di A (A è una matrice quadrata)

• Scrivere lo script sost_indietro.m che implementi il metodo di sostituzione all’indietro per la soluzione di sistemi lineari con matrice dei coefficienti triangolare superiore (eseguire il controllo sulla matrice)

Suggerimento: Indicando con A la matrice dei coefficienti del sistema, porre B=triu(A); e confrontarlo con A

• Scrivere uno script che implementi il metodo di sostituzione

in avanti per la soluzione di sistemi lineari con matrice dei coefficienti triangolare inferiore (eseguire il controllo sulla matrice)

Suggerimento: Indicando con A la matrice dei coefficienti del sistema, porre B=tril(A); e confrontarlo con A

Esercizi

24

Scrivere lo script Matlab esercizio_lab3.m che, assegnata una matrice reale quadrata T: a) calcoli l'elemento di modulo massimo tmax della prima colonna di T; b) se tmax è diverso dal valore assoluto dell'elemento t11, lo porti, con uno scambio di righe, al posto di t11; c) stampi la nuova matrice. Usare come input

Esercizio

25

112

1123

751

T

110

841

32

12

T

40

81

32

T

% Lo script esercizio_lab3 legge in input una matrice quadrata T e scambia la % prima riga con la k-sima riga se si verifica la seguente condizione

% abs(T(k,1)) > abs(T(1,1))

% Se è stato effettuato uno scambio di righe, lo script stampa la matrice T modificata T = input(‘inserire una matrice quadrata, T = ’); [m,n] = size(T); if (m˜=n) error(‘La matrice deve essere quadrata!!!’); end [tmax,ptmax] = max(abs(T(:,1))); if tmax ˜= abs(T(1,1)) % in modo equivalente if ptmax ˜= 1 appo = T(1,:); T(1,:) = T(ptmax,:); T(ptmax,:) = appo; disp(T) end

Soluzione

26

Dal Command Window

>> esercizio_lab3

inserire una matrice quadrata, T = [1 5 7; -3 2 1; 1/2 1 1]

-3.0000 2.0000 1.0000

1.0000 5.0000 7.0000

0.5000 1.0000 1.0000

>> esercizio_lab3

inserire una matrice quadrata, T = [2 1/2 3; 1 4 8; 0 -1 -1]

>>

>> esercizio_lab3

inserire una matrice quadrata, T = [2 3; -1 8; 0 4]

??? Error using ==> esercizio_lab3 at 9

La matrice deve essere quadrata!!!

Soluzione

27

• Scrivere lo script MEG_pivot.m che implementi il metodo di

eliminazione di Gauss con pivoting parziale per la soluzione di sistemi lineari. Usare lo script per risolvere il problema del circuito elettrico

Suggerimento: si modifichi lo script MEG_nopivot.m

adattando lo script esercizio_lab3.m al generico passo k del metodo di eliminazione di Gauss

Esercizio

28

Per usare le istruzioni in esercizio_lab3 al generico passo k del metodo di eliminazione di Gauss (script MEG_nopivot.m) per la soluzione del sistema lineare Ax = b, occorre fare i seguenti cambiamenti, includendo anche lo scambio delle componenti corrispondenti del vettore dei termini noti b [tmax,ptmax] = max(abs(A(k:n,k))); if tmax ˜= abs(A(k,k)) appo = A(k,:); A(k,:) = A(k+ptmax-1,:); A(k+ptmax-1,:) = appo; disp(A) % si può omettere appo = b(k); b(k) = b(k+ptmax-1); b(k+ptmax-1) = appo; end Queste istruzioni vanno inserite nello script MEG_nopivot prima del calcolo del moltiplicatore m (subito dopo il for sull’indice k e prima del for sull’indice i)

Soluzione

29

% meg_pivot % Risolve il sistema lineare Ax=b usando il metodo di % eliminazione di Gauss con pivoting parziale

clear all % cancella tutte le variabili in Workspace %lettura della matrice dei coefficienti e del vettore dei termini noti A = input(‘inserisci la matrice dei coefficienti del sistema A= ’); b = input(‘inserisci il vettore dei termini noti b = ‘’); [m,n] = size(A); % controllo sulla dimensione di A e b if (m˜=n) || ((m==n) && (n˜=length(b))) error(‘Attenzione alla dimensione di A e b!!!’); end % controllo dell’esistenza di un’unica soluzione if det(A) == 0 error(‘la matrice dei coefficienti e’’ singolare!!!’) end

Soluzione (script meg_pivot.m)

30

31

%controllo dei minori principali di testa for i=1:n-1 if det(A(1:i,1:i)) == 0 error('uno dei minori principali di testa e'' singolare!!!') end End for k =1:n-1 % gestisce il numero di passi del metodo di Gauss for i = k+1:n [tmax,ptmax] = max(abs(A(k:n,k))) if tmax ˜= abs(A(k,k))

appo = A(k,:); A(k,:) = A(k+ptmax-1,:); A(k+ptmax-1,:) = appo; disp(A) % si può omettere appo = b(k); b(k) = b(k+ptmax-1); b(k+ptmax-1) = appo; end

m = A(i,k)/A(k,k); % calcolo dei moltiplicatori

32

for j=k:n % modifica degli elementi sulla i-sima riga A(i,j) = A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); % modifica della i-sima riga componente di b end end % metodo di sostituzione all’indietro % All’uscita dei cicli annidati precedenti la matrice A è triangolare superiore x(n) = b(n)/A(n,n); for i =n-1:-1:1 x(i) = (b(i)-sum(A(i,i+1:n).*x(i+1:n)))/A(i,i); end x = x’; % la soluzione deve essere un vettore colonna

fprintf('la soluzione del sistema e'' x= \n ')

fprintf('%+6.5f \n',x(:))

Dal Command Window >>MEG_pivot inserisci la matrice del sistema A= [11 -5 0 0 0 -1; -20 41 -15 0 -6 0; 0 -3 7 -4 0 0; 0 0 -1 2 -1 0; 0 -3 0 -10 28 -15; -2 0 0 0 -15 47] inserisci il vettore dei termini noti b = [500 0 0 0 0 0]' la soluzione del sistema e' x= +70.00000 +52.00000 +40.00000 +31.00000 +22.00000 +10.00000

Soluzione

33

Sistema lineare relativo al problema del circuito senza scambio di equazioni

Dal Command Window >>MEG_pivot inserisci la matrice del sistema A= [11 -5 0 0 0 -1; 0 -3 0 -10 28 -15; -2 0 0 0 -15 47; 0 0 -1 2 -1 0; -20 41 -15 0 -6 0; 0 -3 7 -4 0 0] inserisci il vettore dei termini noti b = [500 0 0 0 0 0]' la soluzione del sistema e' x= +70.00000 +52.00000 +40.00000 +31.00000 +22.00000 +10.00000

Soluzione

34

Sistema lineare relativo al problema del circuito con scambio di equazioni

OSS: >> X=A\b X = 70.0000 52.0000 40.0000 31.0000 22.0000 10.0000

Esercizio • Scrivere uno script Matlab che confronti il tempo richiesto

per risolvere numericamente un sistema lineare usando il solutore di Matlab e il metodo di eliminazione di Gauss con pivoting scritto senza usare il solutore di Matlab

Esercizio

Scrivere uno script Matlab che risolva un sistema lineare con lo algoritmo di Thomas

Algoritmo di Thomas

Algoritmo di Thomas

Algoritmo di Thomas % thomas_diag risolve il sistema lineare Ax=b usando

% l'algoritmo di Thomas, con A matrice tridiagonale

% quadrata avente d come diagonale principale, a e s come

% codiagonali inferiore e superiore

%

% INPUT

% d = diagonale principale di A (vettore)

% a = diagonale inferiore di A (vettore)

% s = diagonale superiore di A (vettore)

% b = termine noto

%

% OUTPUT

% x = vettore soluzione

d = input(‘introduci la diagonale principale di A (vettore) d = ’); a = input(‘introduci la diagonale inferiore di A (vettore) a = ‘); s = input(‘introduci la diagonale superiore di A (vettore) s = ‘); b = input(‘introduci il vettore dei termini noti b = ’);

Algoritmo di Thomas % Controllo delle variabili di input: matrice quadrata

la = length(a);

ls = length(s);

ld = length(d);

if (la~=ls)|(ld~=(la+1))

disp('LA MATRICE NON è QUADRATA!!!')% alternativa al comando error

x = [];

else

n=ld;

tic

% fattorizzazione LU di A

u(1) = d(1);

v = s;

for i = 1:n-1

alfa(i) = a(i)/u(i);

u(i+1) = d(i+1) - alfa(i)*v(i);

end

Algoritmo di Thomas

% calcolo della soluzione L*y = b

y(1) = b(1);

for i=2:n

y(i) = b(i)-alfa(i-1)*y(i-1);

end

% U*x = y

x(n) = y(n)/u(n);

for i = n-1:-1:1

x(i) = (y(i)-v(i)*x(i+1))/u(i);

end

toc

end

whos

Algoritmo di Thomas % thomas.m Risolve il sistema lineare Ax=b con A matrice

% tridiagonale quadrata, usando l'algoritmo di Thomas

% INPUT

% A = matrice tridiagonale

% b = termine noto

%

% OUTPUT

% x = vettore soluzione

A = input(‘inserisci la matrice del sistema A = ’);

b = input(‘inserisci il vettore dei termini noti b = ’);

% Controllo delle variabili di input

% matrice tridiagonale

T = triu(A,-1)-triu(A,2);

B = (A == T);

if ~all(B(:))

error('LA MATRICE NON è TRIDIAGONALE!!!')

end

Algoritmo di Thomas % matrice non singolare

d = det(A);

if d==0

error('LA MATRICE è SINGOLARE!!!')

end

% matrice quadrata

[n,m] = size(A)

if n~=m

error('LA MATRICE DEVE ESSERE QUADRATA')

end

tic

% estraggo la diagonale principale di A e quelle sup e inf

d = diag(A); % diagonale

a = diag(A,-1); % diagonale inferiore

s = diag(A,1); % diagonale superiore

Algoritmo di Thomas % fattorizzazione LU di A

u(1) = d(1);

v = s;

for i = 1:n-1

alfa(i) = a(i)/u(i);

u(i+1) = d(i+1) - alfa(i)*v(i);

end

% calcolo della soluzione

% L*y = b

y(1) = b(1);

for i=2:n

y(i) = b(i)-alfa(i-1)*y(i-1);

end

% U*x = y

x(n) = y(n)/u(n);

for i = n-1:-1:1

x(i) = (y(i)-v(i)*x(i+1))/u(i);

end

toc

Algoritmo di Thomas

Si genera una matrice tridiagonale >> A = diag(ones(N-1,1))+diag(-.5*ones(N-2,1),1)+...

diag(-.5*ones(N-2,1),-1);

>> b = zeros(1,N-1);, b(1) = .5;

>> thomas

inserisci la matrice del sistema A = A

inserisci il vettore dei termini noti b = b

OSS: lo script non ha comandi di visualizzazione del risultato ma la variabile x è nella memoria!

Esercizio • Si consideri il sistema Ax = b, con matrice dei coefficienti A = [4 1 0 0 0; 1 4 1 0 0; 0 1 4 1 0; 0 0 1 4 1; 0 0 0 1 4] e vettore dei termini noti b = [10 12 12 12 10]’. Lo si risolva usando il metodo di eliminazione di Gauss (usando gli operatori predefiniti di Matlab) e l’algoritmo di Thomas e si confrontino i risultati in termini di precisione della soluzione e tempo di esecuzione

Scrivere uno script Matlab che individui gli estremi dell’intervallo di separazione della radice dell’ equazione non lineare f(x)=0 contenuta nell’intervallo I = [a,b] usando opportunamente le funzioni Matlab find e diff. La funzione f e gli estremi a e b devono essere introdotti dall’utente Eseguire lo script per le seguenti funzioni • f() = e + 0.435 (e - 1) / - 1.564 in I=[0.05,0.15] • f(x) = log(x+1)+(x+2)1/2-1 in I=[-0.9,0.5] • f(x) = x 3 – 10x 2 + 5 in I=[0.6,0.8] • f(x) = (x-1)2-3log(x) in I=[0.5,2] • f(x) = sin(x)+x-1 in I=[0,π] • f(x) = log(x)+x2-x in I=[0.7,2.3]

Esercizio

47

% sep_radici_tab % script per la separazione delle radici dell’equazione non lineare % f(x)=0 contenute nell’intervallo I = *a,b] f = input(‘introduci la funzione f = ’) a = input(‘introduci l’’estremo inferiore dell’’intervallo a = ’) b = input(‘introduci l’’estremo superiore dell’’intervallo b = ’) np = input(‘numero di punti in cui suddividere l’’intervallo np = ’) if isempty(np) np=10; end x=linspace(a,b,np); y=f(x); pos = find(abs(diff(sign(y)))); intervallo(1,:) = x(pos); intervallo(2,:) = x(pos+1); fprintf('l''estremo inferiore e’’ %6.5f\nl''estremo superiore e’’ %6.5f\n', intervallo) 48

Esercizio

Dal Command Window >> sep_radici_tab introduci la funzione f = @(x)[x.^3-10*x.^2+5]; f = @(x)[x.^3-10*x.^2+5] introduci l'estremo inferiore dell'intervallo a = 0.6 a = 0.6000 introduci l'estremo superiore dell'intervallo b = 0.8 b = 0.8000 numero di punti in cui suddividere l'intervallo np = np = [] l'estremo inferiore e' 0.73333 l'estremo superiore e' 0.75556 49

Esercizio

Dal Command Window >> sep_radici_tab introduci la funzione f = @(x)[x.^3-10*x.^2+5]; f = @(x)[x.^3-10*x.^2+5] introduci l'estremo inferiore dell'intervallo a = 0.6 a = 0.6000 introduci l'estremo superiore dell'intervallo b = 0.8 b = 0.8000 numero di punti in cui suddividere l'intervallo np = np = 11 l'estremo inferiore e' 0.72000 l'estremo superiore e' 0.74000 50

Esercizio

51

Scrivere uno script matlab che implementi il metodo di bisezione usando un numero fissato di iterazioni Soluzione: versione molto semplificata a = 1; b = 5; for j = 1:10 c = (a+b)/2; if f(c)*f(a) > 0 a = c; else b = c; end end

Esercizio

52

• Si suppone di conoscere la funzione f e che gli estremi a e b siano tali che f(a)f(b)<0

• Il numero di iterazioni è fissato arbitrariamente

• Non si controlla se la approssimazione alla k-sima iterazione è proprio la radice

% bisez_1 % script per la soluzione dell’equazione non lineare f(x)=0 % nell’intervallo I = *a,b] usando il metodo di Bisezione f = input(‘introduci la funzione f = ’) a = input(‘introduci l’’estremo inferiore dell’’intervallo a = ’) b = input(‘introduci l’’estremo superiore dell’’intervallo b = ’) max_it = input(‘numero massimo di iterazioni consentite max_it = ’) for k = 1:max_it c = (a+b)/2; if f(c) == 0 fprintf(‘la radice e’’ %6.5f \n’,c) break end if f(c)*f(a) > 0 a = c; else b = c; end end

Esercizio

53

• la funzione f è richiesta in input ma si suppone che gli estremi a e b siano tali che f(a)f(b)<0

• Il numero di iterazioni è richiesto in input

• Si controlla se la approssimazione alla k-sima iterazione è proprio la radice

Nello script precedente il numero di iterazioni effettuate è il valore che la variabile k assume al termine dell’esecuzione e la approssimazione della radice è contenuto nella variabile c Dal command window >> bisez_1 introduci la funzione f = @(x)[x.^3-10*x.^2+5]; introduci l'estremo inferiore dell'intervallo a = intervallo(1) a = 0.7200 introduci l'estremo superiore dell'intervallo b = intervallo(2) b = 0.7400 numero massimo di iterazioni consentite max_it = 100 la radice e' 0.73460 >> disp(k) 48 >> disp(c) 0.7346

Esercizio

54

Dal command window >> bisez_1 introduci la funzione f = @(x)[x.^3-10*x.^2+5]; introduci l'estremo inferiore dell'intervallo a = .6 a = 0.6000 introduci l'estremo superiore dell'intervallo b = .8 b = 0.8000 numero massimo di iterazioni consentite max_it = 100 la radice e' 0.73460 >> disp(k) 51 >> disp(c) 0.7346 >> disp(f(c)) 0

Esercizio

55

1. Come si può modificare lo script precedente se si vuole introdurre il criterio di arresto a priori?

2. Come si può modificare lo script precedente se si vuole introdurre almeno un criterio di arresto a posteriori?

3. Ricordando le funzioni error e return, modificare lo script precedente usando uno dei due comandi per interrompere l’esecuzione dello script se gli estremi a e b introdotti da tastiera non verificano la condizione f(a)f(b)<0 oppure se uno dei due è proprio la radice cercata. Stampare opportunamente il risultato.

Esercizi

56

1. Come si può modificare lo script precedente se si vuole introdurre il criterio di arresto a priori?

a) max_it va calcolato all’interno dello script nel modo seguente max_it = ceil(log2(b-a)-log2(eps)) b) eps è la tolleranza richiesta e deve essere data come input

quindi la prima parte dello script bisez_1 va modificata come segue f = input(‘introduci la funzione f = ’) a = input(‘introduci l’’estremo inferiore dell’’intervallo a = ’) b = input(‘introduci l’’estremo superiore dell’’intervallo b = ’) eps = input(‘introduci la precisione richiesta sulla soluzione eps = ’) max_it = ceil(log2(b-a)-log2(eps))

Soluzione

57

2. Come si può modificare lo script precedente se si vuole introdurre almeno un criterio di arresto a posteriori? x0 = a; for k = 1:max_it c = (a+b)/2; if f(c) == 0 fprintf(‘la radice e’’ %6.5f \n’,c) break

elseif abs(x0-c)<eps break end x0=c; if f(c)*f(a) > 0 a = c; else b = c; end end

58

Soluzione

E’ necessario definire una variabile in cui si conserva il valore dell’elemento della successione calcolato all’iterazione precedente in modo da poter valutare la differenza tra due approssimazioni successive

3. Ricordando le funzioni error e return, modificare lo script precedente usando uno dei due comandi per interrompere l’esecuzione dello script se gli estremi a e b introdotti da tastiera non verificano la condizione f(a)f(b)<0 oppure se uno dei due è proprio la radice cercata. Stampare opportunamente il risultato. f = input(‘introduci la funzione f = ’) a = input(‘introduci l’’estremo inferiore dell’’intervallo a = ’) b = input(‘introduci l’’estremo superiore dell’’intervallo b = ’) if f(a)f(b)>0 error(‘Attenzione! Non sono verificate le condizioni di applicabilita’’

del metodo di bisezione!!!’) end

59

Soluzione

3. Ricordando le funzioni error e return, modificare lo script precedente usando uno dei due comandi per interrompere l’esecuzione dello script se gli estremi a e b introdotti da tastiera non verificano la condizione f(a)f(b)<0 oppure se uno dei due è proprio la radice cercata. Stampare opportunamente il risultato. if f(a)*f(b)>0 error(‘Attenzione! Non sono verificate le condizioni di applicabilita’’ del

metodo di bisezione!!!’) elseif f(a)*f(b)==0 if f(a)==0 disp(*‘la soluzione e’’ ’, num2str(a)+) else disp(*‘la soluzione e’’ ’, num2str(b)+) end return end

60

Soluzione

61

% bisez_1 % script per la soluzione dell’equazione non lineare f(x)=0 % nell’intervallo I = *a,b] usando il metodo di bisezione f = input('introduci la funzione f = ') a = input('introduci l’’estremo inferiore dell’’intervallo a = ') b = input('introduci l’’estremo superiore dell’’intervallo b = ') if f(a)*f(b)>0 error('Attenzione! Non sono verificate le condizioni di applicabilità del metodo di bisezione!!!') elseif f(a)*f(b)==0 if f(a)==0 disp(['la soluzione e'' ', num2str(a)]) else disp(['la soluzione e'' ', num2str(b)]) end return end eps = input('introduci la precisione richiesta sulla soluzione eps = ') max_it = ceil(log2(b-a)-log2(eps))

62

x0 = a; for k = 1:max_it c = (a+b)/2; if f(c) == 0 fprintf('la radice e'' %6.5f \n',c) break elseif abs(x0-c)<eps break end x0=c; if f(c)*f(a) > 0 a = c; else b = c; end end fprintf('l''approssimazione della radice che soddisfa la precisione richiesta e'' %10.8f \n',c)

63

Esercizio

64

f = @(x)[1./sqrt(x) + 2*log10(0.0001/3.7 + 2.51./(300000*sqrt(x)))] figure, fplot(f,[.0001,.1]) hold on, fplot(@(x)[0.*x],[.0001,.1])

Soluzione: separazione grafica

La radice è contenuta nell’intervallo I =[0.01, 0.02]

Dal command window >> bisez_1 introduci la funzione f = @(x)[1./sqrt(x)+2*log10(0.0001/3.7+2.51./(300000*sqrt(x)))] introduci l’’estremo inferiore dell’’intervallo a = 0.01 introduci l’’estremo superiore dell’’intervallo b = 0.02 introduci la precisione richiesta sulla soluzione eps = 10^(-8) max_it = 20 l'approssimazione della radice che soddisfa la precisione richiesta e' 0.01543061 >>k k = 20

65

Soluzione: script bisez_1

Dal command window >> bisez_1 introduci la funzione f = @(x)[1./sqrt(x)+2*log10(0.0001/3.7+2.51./(300000*sqrt(x)))] introduci l’’estremo inferiore dell’’intervallo a = 0.02 introduci l’’estremo superiore dell’’intervallo b = 0.03 Error using bisezione_script_for2 (line 8) Attenzione! Non sono verificate le condizioni di applicabilità del metodo di bisezione!!!

66

Soluzione: script bisez_1

Se si aggiunge hold on plot(k,c,'b*') nel blocco di istruzioni del ciclo for subito dopo il calcolo di c, aprendo la finestra grafica prima del ciclo, si visualizza anche la successione generata dal metodo. I comandi xlabel, ylabel e title sono stati introdotti subito dopo lo end del ciclo for Dal command window >> bisez_1 introduci la funzione f = @(x)[1./sqrt(x)+2*log10(0.0001/3.7+2.51./(300000*sqrt(x)))] introduci l’’estremo inferiore dell’’intervallo a = 0.01 introduci l’’estremo superiore dell’’intervallo b = 0.02 introduci la precisione richiesta sulla soluzione eps = 10^(-8) max_it = 20 l'approssimazione della radice che soddisfa la precisione richiesta e' 0.01543061

67

Osservazione: script bisez_1

Se si aggiunge un comando di stampa della variabile c ad ogni iterazione, per esempio fprintf(‘%6.5f \n’,c’) subito dopo il calcolo di c nel blocco di istruzioni del ciclo for, si visualizza il seguente risultato 0.01500 0.01750 0.01625 0.01563 0.01531 0.01547 0.01539 0.01543 0.01545 0.01544 0.01543 0.01543 0.01543 0.01543 0.01543 0.01543 0.01543 0.01543 0.01543 0.01543

l'approssimazione della radice che soddisfa la precisione richiesta e' 0.01543061 68

Osservazione: script bisez_1

Cambiare formato di stampa e confrontare i risultati

• Scrivere uno script Matlab che implementi il metodo delle approssimazione successive per la soluzione di un’equazione non lineare usando come criterio di arresto il criterio di arresto a posteriori. Quali parametri devono essere inseriti dall’utente? Stampare il numero di iterazioni eseguite e l’approssimazione prodotta. Stampare e disegnare il vettore contenente la stima dell’errore commesso ad ogni iterazione.

Si considerino le seguenti funzioni di iterazione

relative all’equazione non lineare nell’intervallo • Ripetere l’esercizio precedente usando opportunamente sia il criterio di

arresto a priori che il criterio di arresto a posteriori

• Usare lo script opportunamente per approssimare la radice dell’equazione di Colebrook

Esercizi

69

% script punto_unito % cerca il punto unito xn della funzione f con precisione eps % scegliendo x0 come punto iniziale. Il procedimento iterativo % si interrompe se si raggiunge il numero massimo di iterazioni % maxiter; in caso contrario calcola il numero di iterazioni % effettuate n_iter % % INPUT % f = espressione della funzione di cui si vuole cercare il punto unito % x0 = punto iniziale % eps = limite superiore dell'errore da usare come criterio di arresto % (err<eps con err = |xn-x(n-1)|) % maxiter = numero massimo di iterazioni consentite

Soluzione

70

format long; f = input(‘introduci la funzione di iterazione f = ‘); x0 = input(‘introduci il punto iniziale x0 = ’); eps = input(‘introduci la precisione richiesta eps = ’); maxiter = input('introduci il numero massimo di iterazioni consentite maxiter = '); if maxiter<=0 , error(‘maxiter deve essere un numero positivo’), end figure, hold on, xlabel(‘iterazioni’), ylabel(‘errore’) % grafica errore ad ogni iterazione for n_iter = 1:maxiter xn = f(x0); err = abs(xn-x0); plot(n_iter,err,‘b*') if err<=eps, break, else, x0 = xn; end end fprintf('n_iter \t xn \t errore \n') fprintf('%3d \t %15.10f \t %11.10f \n',n_iter, xn, err) if (n_iter == maxiter) && (err>eps) fprintf('La soluzione non ha raggiunto la precisione richiesta in %5d iterazioni',maxiter)

end 71

Nota: per usare un criterio di arresto a posteriori è stato forzato l’arresto della ripetizione con il comando break

72

>> punto_unito introduci la funzione di iterazione f = @(x)[x.^3+6*x.^2+x-8] introduci il punto iniziale x0 = 1 introduci la precisione richiesta eps = 0.5*10^-5 introduci il numero massimo di iterazioni consentite maxiter = 20 n_iter xn errore 20 NaN NaN

Dal command window

L’errore non tende a zero!!!

73

>> punto_unito introduci la funzione di iterazione f = @(x)[sqrt(8./(x+6))] introduci il punto iniziale x0 = 1 introduci la precisione richiesta eps = 0.5*10^-5 introduci il numero massimo di iterazioni consentite maxiter = 20 n_iter xn errore 5 1.0641779291 0.0000022355

Dal command window

L’errore tende a zero!!! Esercizio: Aumentare la precisione richiesta e graficare l’errore con segno

74

Esercizio

75

Soluzione 1: da Command window

separazione grafica >> g = 9.81; >> u = 2510; >> M0 = 2.8*10^6; >> m=13.3*10^3; >> vi = 335; >> f = @(x)[u*log(M0./(M0-m*x))-g*x-vi]; >> figure, fplot(f,[0,M0/m]) >> hold on, fplot(@(t)[0.*t],[0,M0/m]) >> xlabel('t'), ylabel('f(t)') >> bisez_1 introduci la funzione f = f introduci l’’estremo inferiore dell’’intervallo a = 0; introduci l’’estremo superiore dell’’intervallo b = M0/m-.0001; introduci la precisione richiesta sulla soluzione eps = 0.5*10^-3; max_it = 19

76

Soluzione 1: da Command window

105.26311 52.63155 78.94733 65.78944 72.36839 69.07891 70.72365 71.54602 71.13483 70.92924 70.82645 70.87784 70.90354 70.89069 70.88427 70.88106 70.87945 70.87865 70.87825 l'approssimazione della radice che soddisfa la precisione richiesta e' 70.87824617

77

Soluzione 2: script_razzo.m

close all, clear all % parametri relativi alla funzione di cui si vuole cercare una radice g = 9.81; u = input('inserire la exhaust velocity del razzo u = ') M0 = input('inserire la massa del razzo al momento del lancio M0 = ') m = input('inserire il tasso di consumo di carburante m = ') vi = input('inserire la velocita'' di interesse vi = ') % funzione di cui si vuole cercare una radice f = @(t)[u*log(M0./(M0-m*t))-g*t-vi]; figure, fplot(f,[0,M0/m]) hold on, fplot(@(t)[0.*t],[0,M0/m]) xlabel('t'), ylabel('f(t)') % Richiamo lo script che implementa il metodo di bisezione bisez_1

In questo modo è possibile - studiare il problema al variare dei parametri di input; - confrontare la soluzione prodotta da diversi metodi (inserire le chiamate

allo script punto_unito.m e, dopo averlo scritto, allo script Newton.m

78

Dal Command window >> script_razzo inserire la exhaust velocity del razzo u = 2510 inserire la massa del razzo al momento del lancio M0 = 2.8*10^6 inserire il tasso di consumo di carburante m = 13.3*10^3 inserire la velocita' di interesse vi = 335 inserire l'estremo inferiore dell'intervallo a = 0 inserire l'estremo superiore dell'intervallo b = M0/m-.0001 introduci la funzione f = f introduci l’’estremo inferiore dell’’intervallo a = 0; introduci l’’estremo superiore dell’’intervallo b = M0/m-.0001; introduci la precisione richiesta sulla soluzione eps = 0.5*10^-3; max_it = 19

Soluzione 2: script_razzo.m

79

Soluzione 2: da Command window

105.26311 52.63155 78.94733 65.78944 72.36839 69.07891 70.72365 71.54602 71.13483 70.92924 70.82645 70.87784 70.90354 70.89069 70.88427 70.88106 70.87945 70.87865 70.87825 l'approssimazione della radice che soddisfa la precisione richiesta e' 70.87824617

Funzione print

Per salvare una figura come immagine in uno specifico formato, si può usare il comando

print(fig, nomefigura, formato)

dove fig è la variabile che si riferisce alla figura da salvare

Esempio

>> fig = figure;

>> x = 0:pi/100:2*pi;

>> y = sin(x);

>> plot(x,y)

>> print(fig,'IlMioPlot','-dpng')

80

Salva il grafico come immagine in formato png Il nome del file è IlMioPlot.png

Esercizio

Soluzione

fig1 = figure,

figure, fplot(f,[0,M0/m])

hold on, fplot(@(t)[0.*t],[0,M0/m])

xlabel('t'), ylabel('f(t)')

print(fig1,‘separazioneGrafica','-dpng')

81

Modificare lo script precedente salvando la figura relativa alla separazione grafica in formato .png

Oss: Il file separazioneGrafica.png viene salvato nella directory di lavoro