Appunti Metodi numerici per l'ingegneria navale

Esame

Orale, da portare il PC perché ci viene chiesto di fare un esercizio su Matlab e ci fa domande di teoria e controlla esercitazioni (solo codice).

Comandi di base

• Clear: cancella workspace.
• Salvare nome senza spazi e iniziare con lettere!
• Close all: cancella finestre grafiche.
• Se si mette ";" non compaiono i risultati nel workspace.

Vettori e matrici

Vettori e matrici

V = [1 4 3 2] vettore riga.

Posso usare spazi o virgole per separare elementi in una riga (in Python uso le virgole).

LV = V' Trasposta (in questo caso, vettore colonna) - SV = [1;4;5;7] Uso ; per dividere elementi di una colonna. Definisco un vettore con [] 134.

A = [1,3,4; 5,7,8;1,1,2] > Richiamo elemento con () 578112.

INDICEDOLORE~ [2]

a = V(3) ↪ Lettore 5.

a = 134.

a = A(1,2) A 3.

a => =578=LCONN 112.

Sottomatrici

Voglio prendere i valori delle righe da 1 a 2.

Sottomatrice B = A(1:2, 1:2) 13413A B-518 57=.

Voglio prendere i valori delle colonne da 1 a 2 =112134 13B = (:, 1:2) A B>578= = = 57↪ 112.

Righe 11 Tutte L Db = B (end ; 2) B 1==Lultimo elemento 134 =.

Voglio prendere i valori di tutte le righe.

B = A (:, [1, 3]) A B578 =>=.

Voglio prendere i valori della 1^{^} e 3^{^} colonna 112.

Help

Se scrivo nella workspace "doc nomecomando" si apre online l'help di un comando.

Operazioni matrici

Alcune funzioni restituiscono un vettore di dati, es. Size restituisce n° righe e n° colonne di una matrice.

Posso utilizzare queste formattazione.

N (B) Size= -noriGHe1 [3 2) NB = == ; no colonne ↪12 (B)b)[a b 23Size a ===: % =Fzeros(3 2) eve(3)C F= =;= = =-888 [A B]=3)(3 GD ones ; === ,000 ↓=337.

Con la "," le affianco mentre * 2)a 3 Ones(2 EE ,a =>== ";" le incolonno [devono concordare righe e colonne].

• Command window: dove posso scrivere comandi direttamente.
• Workspace: mostra le variabili attive.

Grafica su Matlab

Non sono da sapere i Close all chiude tutte le finestre. Clear comandi a memoria, Clear cancella workspace. Close all usare "help". Clc cancella command window.

X = linspace(0, 2*pi, 150);

Y1 = sen(x); x vettore con un certo numero di punti equispaziato.

Y2 = cos(x); da 0 a 2pigreco voglio 150 punti.

Plot(x,y1,'r';x,y2'g') Per aver la stessa scala in x e y.

Legend('sin','cos')

Axis equal gcf è la figura corrente (salva tutto della figura).

Saveas(gcf,'figura1.png') gca sono gli assi correnti (salva l'interno).

Grid on Grid on abilita la griglia.

Secondo modo: Figure: apre una nuova finestra grafica.

Figure Una prima modifica è fare plot(x,y,'r:').

Plot(x,y1,'r') ovvero metto la curva rossa e a puntini.

Hold on.

Plot(x,y2, 'g') Se metto ; a fine della riga non mi compare.

Legend('sen', 'cos') il comando nella command window.

Axis equal.

Grid on.

Copio il primo grafico e scrivo hold on, questo continua a plottare.

Hold off sulla stessa finestra. Rimane attivo fino a che non metto hold off, dopo questo ogni nuovo grafico sostituirà quello precedente (la finestra rimane aperta ma il contenuto viene sovrascritto).

Se volessi aprire un'altra figura e non sovrascrivere devo mettere di nuovo "figure".

Esempio da help "plot"

x = -pi:pi/10:pi;

y = tan(sin(x)) - sin(tan(x));

tratteggiatalinea, figure -T monkerGreen cone "..." servono per mandare a capo il comando, si fa per mettere le Square.

"plot(x,y,'--gs',... cose ordinate (spessore 'LineWidth',2,... linea) 'MarkerSize',10,... 'MarkerEdgeColor','b',... 'MarkerFaceColor',[0.5,0.5,0.5]) del marker Interne~ La tripletta RGB.

Funzioni

FUNZIONI

Perché bisogna fare sottoprogrammi? (Funzioni)

• Uso frequente di una stessa parte di codice.
• Mantenimento e lettura del codice.
• Nome delle variabili.

Iniziano con la parole chiave function. Vanno messe in fondo al file o in un file separato. Ad esempio:

Function scrivi(l) l è una stringa di caratteri.

Disp(l)

end

Le function vanno chiuse con end.

Per eseguirla va copiata tutta la prima riga ma senza "function" e lo mettiamo nel programma principale.

l='ciao';

scrivi(l)

Altro esempio

x=3;

a=2;

Function output = nome(input)

y= f1(x,a) L'output y deve avere un valore alla fine.

Function y= f1(x,a)

y= a*x²;

End

x=3;

a=2; Funziona lo stesso perché conta la posizione dei parametri di input.

b=6; e output e non il nome!

y= f1(x,b)

Le variabili in una funzione sono locali (in qualunque linguaggio di programmazione).

Function y= f1(x,a)

y= a*x²;

End

Sistema lineare

Creare una funzione che prende in input una matrice A e un vettore colonna b e restituisce la soluzione del sistema A*x=b.

Se mettessi la funzione in un A = [1,2;6,3];

A = [1,2; 6,3] b = [3;2];

altri file devo chiamare il file.

b = [3; 2] [x,c] = sistema(A,b).

x = sistema(A,b) con il nome della funzione.

Function [x,c] = sistema(A,b)

Function x = sistema(A,b) x = A^-1*b;

Per avere due uscite:

x = A^-1*b · c=2;

End (vettore in uscita) End.

Lunghezza, minimo e massimo

Creare una funzione che prende in input un vettore v e restituisce in output la lunghezza del vettore, il minimo e il valore massimo.

v = [1,4,7,-1,2.3]; non è importante se la chiamo sopra in maniera.

[lunghezza,minimo, massimo]= casa(v) diversa rispetto a come l'ho chiamata nella funzione.

Function [l,m,M]= casa(v); Se usassi una matrice questo mi rida la riga.

l=length(v); massima o minima quindi per avere il singolo.

m=min(v); valore devo usare min(min(v)) o max(max(v)).

M=max(v);

End

Strutture di controllo

Istruzione condizionale if

Esempio: risolvi equazione di secondo grado: a*x² +b*x +c =0.

Ogni volta che metto if devo chiudere con end.

delta = b² -4a*c ;

If delta >0 Solo else è per dire altrimenti, ovvero

disp('due soluzioni reali e distinte') se nessuno dei casi precedenti.

Elseif delta<0

disp('due soluzioni complesse') Elseif e else sono facoltativi.

Elseif delta ==0

disp('due soluzioni reali coincidenti') Xkè tortierenelleda c'èinternetè cercare su non.

End > >= < <= ~= (sarebbe ≠).

X1= (-b-sqrt(delta))/(2*a)

X2= (-b+sqrt(delta))/(2*a) == è per dare la condizione mentre =0 è l'assegnazione.

nota: Matlab conosce il numero complesso i e j.

Cicli

Ci sono due tipi di cicli: ciclo for e ciclo while. Si usa for quando si conosce il numero di iterazioni. Si usa while quando non so conosce il numero di iterazioni e si ferma il ciclo con un criterio di arresto.

Vogliamo sommare i numeri da 1 a N.

N=100; Faccio un vettore da 1 a 100.

V= 1:100; [ ] in disp serve a concatenare due stringhe.

a= sum(V);

disp(['somma =', num2str(s)]) ↪ trasforma un numero in una stringa di caratteri.

Faccio una funzione somma:

s = somma (N)

disp(['somma =', num2str(s)])

Function s=somma(N); Assegno valore iniziale 0 e poi creo un ciclo.

s= 0

For n=1;N Esempio: For n =[1 4 7]

s= s+n A destra dell'= devo avere un vettore.

End

End

Somma dei numeri dispari

Per sommare solo numeri dispari: in Matlab è consigliato usare il più possibile i vettori.

s = somma (N) I cicli rallentano l'esecuzione.

disp(['somma =', num2str(s)])

Function s=somma(N); Mettere tic toc consente di calcolare il tempo (tra tic o toc).

s= 0

For n=1:2:N 2x 1da tic.

Che vas= s+n, diSaltando 2, s = somma (N).

Fino Na toc.

End da 0partopari.

Per i disp(['somma =', num2str(s)])

End

FizzBuzz

Implementare un codice che per tutti i numeri da 1 a N, fa il display a video del numero stesso, ma se è divisibile per 3 fa il disp di "fizz", se è divisibile per 5 fa il disp di "buzz", se è divisibile per 3 e per 5 fa il disp di "fizzbuzz". Usare la funzione mod per vedere se è divisibile per un numero intero.

N = 50; Se ho due condizioni posso mettere and o or.

for n = 1:N

if mod(n, 3) == 0 && mod(n, 5) == 0 In Matlab and si scrive & mentre or |.

disp('FizzBuzz')

elseif mod(n, 5) == 0 And: condizione globale è vera se entrambe.

disp('Buzz') sono vere.

elseif mod(n, 3) == 0 Or: condizione globale è vera se una delle due.

disp('Fizz') else è vera.

disp(num2str(n))

end

end

Fare una function che prende in entrata n e restituisce in uscita la stringa corrispondente all'esercizio precedente: stringa deve essere il numero stesso o fizz o buzz o fizzbuzz.

Stringa = giocoFizzBuzz(n)

Come scorrere elementi di un vettore?

V= [1,5,6,2,15,8,3]; For i=1: length(v)

For i=1: length(v) Va da 1 alla lunghezza del vettore.

Stringa= giocoFizzBuzz(v(i))

Disp(stringa) V Oppure n=v(i) e poi fare stringa=giocoFizzBuzz(n).

End elemento nella posizione i del vettore v.

Non conta il nome della variabile ma solo l'ordine!

Ciclo while

N = input('Dammi N positivo: ');

while N <= 0 Continua fino a quando quella condizione non è vera.

N = input('Dammi N positivo: ');

end Funziona or perché io voglio uscire dal ciclo.

N = input('Dammi un numero positivo e divisibile per 3: ');

while N <= 0 | mod(N, 3) ~= 0 solo quando entrambe.

N = input('Dammi un numero positivo e divisibile per 3: '); le condizioni siano false.

end

disp(['Il numero inserito è: ', num2str(N)]);

While non ha un contatore. Nel caso lo devo creare, es:

i=0

N = input('Dammi un numero positivo e divisibile per 3: ');

while N <= 0 | mod(N, 3) ~= 0

i=i+1

N = input('Dammi un numero positivo e divisibile per 3: ');

end

disp(['hai sbagliato', num2str(i), 'volte']);

Approssimazione di funzioni

Immaginiamo di avere dei dati sperimentali (file "aletta_sez_rettangolare") e voglio sapere il risultato corrispondente per dei valori che non ho.

Nell’ingegneria si ha spesso un modello matematico predefinito che dipende da parametri che non sono noti e che devono essere determinati sperimentalmente. Altre volte, invece, i valori dei parametri del modello sono disponibili, ma devono essere confrontati con i dati sperimentali per verificare la validità del modello stesso o testare la validità di diversi modelli alternativi per stabilire quale sia il migliore in base alle osservazioni empiriche.

Esistono diversi metodi per l’approssimazione di funzioni.

lineare a tratti.

INTERPOLAZIONE polinomiale.

• Approssimazione lineare a tratti: i valori sono approssimati tramite delle rette che congiungono tali valori; questo è il modo con cui ad esempio Matlab crea i grafici di funzioni unidimensionali.
• Interpolazione polinomiale: dati n + 1 valori (xj, yj) con j = 0, . . . , n, esiste un unico polinomio di grado n passante per tali punti. Tale polinomio può essere determinato nel modo più semplice imponendo che il polinomio generico p(x) = a + a x + a x +. . . + a x passi per i valori dati. All’aumentare del numero di 122"no punti, però, può accadere che diventi instabile (fenomeno di Runge). Es7 Pr(X)! punti grado i passante per n.

Quando usare l'interpolazione

Y ...yX · ·1.

1. Quando ho tanti dati.
2. Quando gli errori sono piccoli.

Yo ~Xo :5 2. 81X1· : ho vuote. Se aone potrebbe oscillare AnXn punti)(zone pochi con X eb·X Τ ....~I.

Su Matlab

• interp1 (1D)
• polyfit
• interp2 (2D)

Per leggere un file Excel:

A=readmatrix('Aletta_Sez_Rettangolare.xlsx', 'sheet', 2)

X= A(:,1); % leggo la temperatura.

Y=A(:,2); % conducibilità.

xt= 320; % temperatura di cui voglio conoscere la corrispondente conducibilità.

yt = interp1(X,Y,xt); Il file deve essere nella cartella di lavoro altrimenti.

plot(X,Y,'-or') linea continua Matlab non lo legge. Potrebbe essere anche che non lo,

hold on ↪ circoletti rosso,

plot(xt,yt,'*b') legga se il foglio Excel è aperto. Per leggere il secondo.

hold off foglio di lavoro con le tabelle dobbiamo indicarlo.

Regressione lineare e minimi quadrati

Supponiamo i dati siano più o meno su una retta, voglio trovare la retta che approssima il meglio possibile i punti (non che passa per i punti). "regressione lineare" obaxy +=.. "retta di regressione" vuol dire calcolare a e b.

9319 Impongo che passi per i punti dati: ·· S baXo + Yo= SISTEMAXo 1 Yobax1 FORMAINy1 %+ aX1= SOURADETERMINATO·%· bMATRICIALE poche incognite e: InXn 1 Iequaziontantebaxn yn+ = A BX#Ax risolvibileèB non= VAX-B 0=-m 1n +RvettoreE.

Se non posso fare in modo che Ax-B sia = 0, posso fare in modo che Ax-B sia il più piccolo possibile.

METODI MINIMI QUADRATI Matlab:2 211011Ax-Bmin comando=> [LEAST SQUARES] (regressione) polyfiti 1=-Rn misuro la norma(a b) ↓, Ottimizzazione e eucidea evitare la.

Per edapplicaradice funzione Cla è radicequadratoalnorma una può cambiare) ilquindi mincrescente non.

Ottimizzare vuol dire calcolare il minimo o il massimo di un funzionale: F(X(minXERn.

Condizione necessaria non sufficiente: F(x)=0.

BI(2)(Il AX BAX0 >=- = ATBATAX = 1- pseudoinversaATB(ATA)X =.

Immaginiamo ora i dati su circa una parabola zacalcolare a,b,c. Imposto di nuovo che passi per i punti dati X bxo+ ↑+Coa YoQD :=A BX"ATBla soluzione è sempre CATA)X =.

Si ua il metodo dei minimi quadrati quando:

• polinomio grado n non conoscete la funzione > + calcolare i parametri;
• conoscete la funzione asen(x) bos(x) y += Scartiminaxi yminllAX-BI erb-+= quadraticimedi;
• La b-ysistema retta con + La resta perpassa, punti ili megliopossibile"j.

Esercizio dei minimi quadrati

Per l'esercizio dopo (es. dei minimi quadrati): Voglio calcolare a e b a.

Sappiamo che per la densità abbiamo una funzione del tipo bxy = ai minimi quadrati.

Facciamo come prima e imponiamo il passaggio per i punti dati non posso farlo!

nonuneri è non variante a è l'esponente nei parametri.

Non è un sistema lineare.

minlly-bxa.

Posso ancora dire che minimizzo: R2 (sono (a b) incognite le mie sempre e, minlly-sidu.

Più in generale posso ancora dire che se ho una funzione: f(X c)bay >= ,, , R3(a c)b =, , ~.

La differenza rispetto a prima è che non ho più la è di nuovo un problema di ottimizzazione di un funzionale min F(x)X formula della pseudoinversa (formula analitica), questo problema va risolto numericamente.

bx9.

Linearizzo: y = log(bxP)

log(y) = +b)(4) +a) ret=N YY.

1. Applicare log per calcolare xtilde, ytilde.
2. caloclare a,b ai minimi quadrati con polyfit.
3. a=a; b=exp(b)
4. plot di y =bx (no polyval perché valuta un polinomio e questo non lo è).

Codice su Matlab

grado del polinomio.

n=2;

xt = 320; il potevo mettere direttamente esstringax='temperatura'; perché 'conducibilitàtemperatura' nome, funzione.

L'ho ciscritto nella parolastache unastringa.

y='conducibilità'; Mateab comunque.

Quindi capisce.

[p,yt]=polinomio(X,Y,n,xt,stringax, stringay); L ali ilcalcolarevoglio valorepunto in.

function [p,yt]= esponente(X,Y,xt)

Xtilde = log(X); % per il calore specifico invece: terza colonna. Bisogna.

Dtilde = log(Y); usare un polinomio di grado 2. Basta cambiare n e la.

Ptilde = polyfit(Xtilde,Dtilde); colonna Y=A(:,3) calore specifico.

atilde = p(1);

btilde = p(2);

logD = atilde*Xtilde+btilde; % densità: quinta colonna. Funzione y=b*x^a (da teoria).

a = atilde;

b = exp(btilde); polyfit calcola un polinomio ai minimi quadrati.

y = b * X^a; Perché interpolazione polinomiale e minimi quadrati si fanno.

figure con polyfit? La differenza sta nel fatto che ad esempio se ho.

plot(X,y, 'or'); 10 punti e scelgo grado 9 allora fa l'interpolazione.

hold on

plot(xt,yt,'b*') polinomiale, se invece scegliessi gradi inferiori allora fa il.

hold off minimi quadrati (se non ha abbastanza punti).

end

function [p,yt]=polinomio(X,Y,n,xt,stringax, stringay)

p=polyfit(X,Y,n);

a=p(1); % primo membro vettore p.

b=p(2); % secondo membro vettore p.

x1 = linspace(X(1), X(end), 150);

% y1=a*x1+b;

y1=polyval(p,x1);

figure

plot(X,Y, 'or',x1,y1,'b')

yt = polyval(p,xt); % ovvero 'prendi il polinomio e valuto nel punto xt'

hold on

plot(xt,yt,'b*') polyfit se non ha abbastanza punti fa i minimi quadrati. se

xlabel(stringax)

ylabel(stringay) abbiamo n+1 punti prende un pol