Wirgilio

Scrivo questo breve articolo per suggerirvi qualche trucchetto che può semplificarvi la vita lavorando con MATLAB.

Cerca su Google da command window

Molto spesso mi capita di dover cercare rapidamente qualcosa su Google, mentre lavoro nella console di MATLAB. Specialmente quando lavoro sul netbook, dove ho davvero poca RAM disponibile, aprire Chrome per una rapida ricerca è davvero un supplizio.
Allora ho pensato: perché non utilizzare il browser integrato in MATLAB?
Detto fatto.

Ecco una semplice funzione che richiama il browser integrato per effettuare ricerche immediate:
Download google.zip
Mostra codice ▼

function google(search, varargin)
     % Cerca qualcosa su google nel browser di matlab
     % google('wirgilio')
     % La lingua predefinita è l'inglese, è possibile scegliere una lingua
     % diversa inserendola come secondo argomento
     % google('wirgilio','it')
     if nargin == 1
         lang = 'en';
     else
         lang = varargin{1};
     end
     web(sprintf('http://www.google.com/search?hl=%s&q=%s', ...
     lang, urlencode(search)))
 end
 

In questo modo basterà digitare google(‘oggetto da cercare’) per aprire il browser di MATLAB.

.

Avvia il debug durante l’esecuzione

Durante l’esecuzione di script molto lunghi è fondamentale essere sicuri che il codice stia funzionando correttamente. MATLAB mette a disposizione un comando molto utile che permette di investigare durante l’esecuzione del codice: keyboard.

Il comando è molto semplice è intuitivo, ma il problema è: come faccio ad attivarlo quando decido io? Ad esempio, sto osservando i risultati di uno script ogni tot iterazioni attraverso una figura, quando mi accorgo che qualcosa sta andando storto e voglio avviare il debug. Come fare?
L’idea è quella di utilizzare una variabile presente all’interno delle figure, ossia CurrentCharacter. Quando si preme un tasto della tastiera attraverso una figure, la variable CurrentCharacter assume il valore pari al tasto premuto.

Con una semplice riga si può quindi avviare il debug premendo un tasto della tastiera, in questo caso la lettera k.

if get(gcf,'CurrentCharacter') == 'k', keyboard, end

Per continuare l’esecuzione sarà sufficiente chiudere la figura e premere F5.

Esempio di utilizzo:
Mostra codice ▼

 for i = 1:100
    % ... Calcoli complicati ...
     plot(sin((0:.1:2*pi )+ i))
     title('Premi k per il debug!')
     % Riga di debug
     if get(gcf,'CurrentCharacter') == 'k', keyboard, end
     % % %
    pause(0.1)
 end

.

Focus delle figure

Quando una figure di MATLAB viene selezionata, questa acquisice il focus e viene mostrata davanti a tutte le altre finestre. Se si sta lavorando al codice durante l’esecuzione, può succedere che il continuo apparire delle figure ci impedisca di fare qualsiasi altra cosa.

Una possibile soluzione è quella di utilizzare la funzione set per cambiare la figura attiva, invece del classico figure().

Esempio di utilizzo:
Mostra codice ▼

 f1 = figure(1);
 f2 = figure(2);
 for i = 1:100
     % ... Calcoli complicati ...
     set(0, 'CurrentFigure', f1)
     plot(sin((0:.1:2*pi )+i), 'r')
     title('Figura 1')

     set(0, 'CurrentFigure', f2)
     plot(sin((0:.1:2*pi )-i), 'b')
     title('Figura 2')

     pause(0.1)
 end

.

Figure a schermo intero

Un altro comando che utilizzo molto spesso è quello che permette di aprire una figure a schermo intero. Il trucco è definire le unità normalizzate quando si crea la figura:

figure('Units','normalized','Position',[0 0 1 1])

L’attributo Position richiede un array di quattro valori: [Posizione_x, Posizione_y, Larghezza, Altezza]. Utilizzando le unità normalizzate sarà sufficiente inserire valori unitari per larghezza e altezza.
Questo tipo di unità permette anche una più agevole disposizione delle figure sullo schermo, indipendentemente dalle sue dimensioni.

Esempio di utilizzo:
Mostra codice ▼

f1 = figure('Units','normalized','Position',[0 0 .5 .5]); imagesc(magic(10))
 f2 = figure('Units','normalized','Position',[.5 0 .5 .5]); spy;
 f3 = figure('Units','normalized','Position',[0 .5 .5 .5]); contourf(gallery(3))
 f4 = figure('Units','normalized','Position',[.5 .5 .5 .5]); plot(rand(1,10))
 pause
 close all

.

Barra di caricamento

Una delle funzioni che adoro di MATLAB è senza dubbio la waitbar.
Ideale quando un codice richiede più di qualche istante per essere eseguito!

L’utilizzo è estremamente semplice, riporto direttamente l’esempio dall’help di MATLAB:

h = waitbar(0,'Please wait...');
for i=1:1000,
% computation here %
waitbar(i/1000,h)
end

Se le iterazioni sono molto rapide, come nel codice dell’esempio, la chiamata a waitbar può essere più lenta del calcolo che viene effettuato durante l’iterazione stessa e può rallentare drasticamente l’esecuzione del codice! Per evitare ciò, vi consiglio di utilizzare sempre una condizione che aggiorni la waitbar un numero sensato di volte:

h = waitbar(0,'Please wait...');
 for i=1:1000,
     % computation here %
     if mod(i, 100) == 0
        % Questo codice viene eseguito solo ogni 100 iterazioni
        waitbar(i/1000,h)
     end
 end

Condividi:
Facebook
Twitter
E-mail
Stampa
								

Programmare è un’attività davvero molto stimolante. È un po’ come risolvere i cruciverba o i giochi di logica della Settimana Enigmistica, solo che hai molta più libertà e hai a disposizione strumenti talmente potenti da avere solo la fantasia come unica limitazione.

Non si può dire di saper  usare un computer, se non si conosce un linguaggio di programmazione!

È proprio da questi presupposti (abbinati all’esaurimento-esami) che nasce questo semplice giochino, realizzato con MATLAB, Flood.
L’idea non è mia, ma non credo si possa parlare di un vero e proprio ideatore di questo gioco… È un po’  come giocare con il Paint!

Il gioco è molto semplice: bisogna riempire l’intera finestra di gioco con uno stesso simbolo, nel minor numero di mosse possibile, partendo dall’angolo in alto a sinistra. I simboli simili si contagiano orizzontalmente e verticalmente. È possibile modificare i simboli presenti nel gioco e il livello di difficoltà.
Una matrice generata casualmente, color, contiene i numeri a cui sono associati i simboli, e un vettore di nome path contiene tutti gli elementi che sono stati “riempiti”.

Il gioco non ha un’interfaccia grafica, per il semplice motivo che l’ho realizzato sul tablet con Octave, che ha un po’ di problemi con gli oggetti figure. Sarebbe carino realizzarne una versione grafica utilizzando una matrice di colori, ma è un gioco talmente stupido che non credo ne valga la pena! Ad ogni modo sarebbe sufficiente utilizzare la funzione imagesc applicata alla matrice color, con un’opportuna colormap.
Probabilmente è stato più divertente realizzarlo che giocarci successivamente 🙂

Download – flood.zip (pochi kB)

Download – flood.zip (pochi kB)

Premessa: la lettura di questo articolo potrebbe risultare molto confusa. È talmente specifico che, molto probabilmente, servirà a una persona su un milione. In realtà servirà molto più a me per fissare certe cose che a qualcun altro…

L’induzione magnetica (e la velocità indotta da un vortice, in aerodinamica) è descritta dalla legge di Biot-Savart:

Dove k è una costante moltiplicativa proporzionale all’intensità (della corrente o del vortice). Per questa legge esistono diverse soluzioni analitiche, per i casi più semplici e anche più frequenti: filo rettilineo di lunghezza infinita, solenoide, …

Potrebbe essere interessante osservare il campo magnetico anche per un filo di una forma più strana, come una spirale a passo non costante, o una forma geometrica complessa. Soluzioni analitiche per forme geometriche complesse non ne esistono ed è necessario integrare numericamente la legge di Biot-Savart.
Concettualmente è molto semplice: divido il filo in piccoli intervalli finiti di lunghezza , per ognuno dei quali calcolo l’induzione e infine ne faccio una sommatoria:

Ecco cosa fa il codice che ho scritto: si assegna un dominio di calcolo (non troppo grande, altrimenti ci impiega un’eternità!) e si fa un ciclo per ogni punto del dominio (tre cicli for su i, j, k). Per ognuno di questi punti, si discretizza la curva in intervalli di ampiezza ds e si calcola l’induzione di ogni elementino sul punto del dominio in questione, sommando i vari contributi (componenti U, V, W).

Non è il massimo dal punto di vista dalle performance, è stata la prima cosa che mi è venuta in mente e non ho neanche pensato se fosse la migliore. Funziona e basta 😀

Un unico appunto lo vorrei fare sull’utilizzo delle meshgrid e del prodotto vettoriale. Il campo magnetico (campo di moto) è descritto da una serie di vettori tridimensionali. Questi vettori dipendono dal punto di applicazione, o meglio, ogni punto del campo ha un suo vettore di induzione, che ne descrive direzione, intensità e verso. In MATLAB, ogni cosa è una matrice e non fanno eccezione i campi vettoriali. Il problema è questo: ogni vettore ha tre componenti, ognuna delle quali è contenuta in una matrice NxNxN. Ho bisogno di tre matrici NxNxN per descrivere il campo vettoriale. Come faccio però a capire in quali punti sono applicati questi vettori?
Qui interviene la funzione meshgrid, che genera una matrice contenente le coordinate dei punti in cui voglio che siano applicati i vettori. Quindi avrò: tre matrici U, V e W contenenti le componenti lungo gli assi x, y e z dell’induzione e tre matrici X, Y e Z contenenti i punti di applicazione dei vettori. Quando si fa il calcolo del di Biot-Savart bisogna tenere presente quanto detto: per il calcolo del vettore R dovrò considerare le matrici X, Y e Z, dato che mi interessa il punto di applicazione del vettore.

Forse era una precisazione banale. Forse no.

Considerazioni aerodinamiche

Il codice permette di integrare la legge di Biot-Savart, che è nota nella fisica per l’induzione di un campo magnetico. In aerodinamica viene invece utilizzata per calcolare la velocità indotta da un vortice. È interessante andare a calcolare il campo di moto nella scia di un’ala utilizzando questo codice.

La Teoria del Filetto Portante di Prandtl ci spiega come sia possibile modellare il campo di moto attorno ad un’ala mediante un sistema di vortici. La schematizzazione più semplice di questo sistema di vortici è il cosiddetto vortice a staffa, di cui rubo senza pietà un’immagine da Google:

Inserendo nel codice i punti di un vortice a staffa:

L = [-1 6 0;
     -1 0 0;
      1 0 0;
      1 6 0];

Si ottiene un campo di moto di questo tipo:

Non è molto preciso perché è stato calcolato su un dominio di 20 punti, ma rende comunque l’idea.
Se a questo campo di moto sommiamo una corrente uniforme, otteniamo il campo di moto attorno a un’ala:

Il confronto con un’immagine ottenuta per via sperimentale (linee di fumo) mostra la validità di questa teoria:

Download dell’animazione in 3D (.avi, 4.9 MB)

Un’altra figura interessante che è possibile ottenere è quella relativa al downwash:

L’immagine è stata ottenuta utilizzando la funzione surf alla componente verticale della velocità.

Da questa figura si possono notare due cose: la prima è che il vortice parallelo all’ala (vortice aderente) produce un upwash sulla corrente incidente, la seconda è che il downwash è dovuto essenzialmente all’effetto dei vortici laterali (vortici liberi).

Tutte le figure delle linee di corrente sono state ottenute utilizzando la funzione streamtube, che prende in input le mesh, il campo di moto e tre vettori/matrici contenenti i punti da cui devono partire le linee di corrente.

Potete scaricare il codice dal seguente link, o dal file exchange di MATLAB.

Download biot_savart.zip
Download da MATLAB file exchange

Stavo facendo una simulazione in CFD e volevo ottenere un video a partire da un migliaio di immagini. In giro ci sono decine di software per farlo, ma a prima vista non ho trovato nulla di estremamente rapido, gratuito e che non richiedesse installazione. Allora ho pensato: “Io ho MATLAB!“. Bene, perché non usarlo per questo scopo? Detto fatto. Dopo un’oretta di codice per capire come utilizzare la classe VideoWriter mi è uscito questo simpatico codice che permette di:

1. Selezionare le immagini (in qualsiasi formato)
2. Selezionare un file di output .avi
3. Impostare il framerate, la qualità (JPEG) e un fattore di scala per ridimensionare le immagini
4. Esportare il video

Il tutto con una simpatica interfaccia grafica 😀

*** Download time-lapse.zip ***

Ho scritto il codice in inglese per poterlo caricare sul file exchange di MATLAB.

*** Download time-lapse.zip ***

Giusto per curiosità, ecco il video che ho ottenuto con lo script: http://www.youtube.com/watch?v=xYUQxIDwDOQ

– Download del codice – eulero.zip

Ultimamente ho passato un po’ di tempo a giocare con la funzione ode45 di MATLAB e, dopo aver speso una decina di minuti necessari a capire come impostare i vari parametri di funzionamento, ho iniziato ad integrare tutte le equazioni differenziali che ho incontrato nel corso degli studi.

Tra i vari problemi in cui mi sono cimentato ultimamente, c’è quello delle Equazioni di Eulero, ossia le equazioni che descrivono le rotazioni di un oggetto di forma arbitraria soggetto a dei momenti forzanti. In questo articolo voglio proporre un codice scritto da me che permette la risoluzione di queste equazioni e ne mostra i risultati sottoforma di animazione.

Impostazione del problema

Per prima cosa, impostiamo il problema. Abbiamo un corpo di forma arbitraria nell’ipotesi semplificativa, grazie a opportune simmetrie, di prodotti d’inerzia nulli; chiameremo i momenti principali d’inerzia A, B e C. Questo corpo ha un orientamento nello spazio che è descritto da una terna di angoli di Eulero; in particolare, è stata scelta una sequenza del tipo XYZ e gli angoli sono stati rinominati. Il corpo è soggetto ad un momento le cui componenti attorno agli assi x, y, z, sono state rinominate L(t), M(t), N(t). Le equazioni da risolvere, in questi termini, sono:

La funzione che utilizziamo per integrare le equazioni è, come preannunciato, ode45. Questa funzione accetta in input una function(t,y), il cui valore è il secondo membro di un’equazione differenziale del primo ordine. Il nostro sistema è del secondo ordine, ma è facilmente trasformabile in un sistema del primo ordine utilizzando altre tre equazioni. Chiamando:

Il sistema diventa:

Il codice

– Download del codice – eulero.zip

Per quanto riguarda il codice, questo è diviso in tre file:

• main.m – contenente il codice principale del programma
  
• eulero.m – contenente il sistema di equazioni da integrare
• ang2dcm.m – funzione che restituisce la matrice dei coseni direttori utilizzata per le rappresentazioni grafiche

Iniziamo dal main.m.

Mostra codice ▼

clc; clear; close all
%% Equazioni di Eulero
tfin = 5; % Tempo di simulazione
prec = 1e-4; % Precisione integrazione

deg2rad = 180/pi;

% Scegli il moto
user = menu('Scegli un caso:','Piccolo disturbo','Moto retrogrado','Not a Disk','Spirale');
switch user
 case 1
 A = 1; B = 1; C = 30; % I_x, Y_y, I_z
 L = @(t) interp1([0 0.2 0.3 0.4 tfin], ...
 [0 0 5e3 0 0 ], t, 'linear');
 M = @(t) 0;
 N = @(t) 0;
 INIT = [0; 0; 50; 0; 0; 0]; % omega_x, omega_y, omega_z, psi, theta, phi
 case 2
 A = 1; B = 1; C = 10;
 L = @(t) interp1([0 0.2 0.3 0.4 tfin], ...
 [0 0 7e2 0 0 ], t, 'linear');
 M = @(t) 0;
 N = @(t) 0;
 INIT = [0; 0; 5; 0; 0; 0];
 case 3
 A = 1; B = 5; C = 10;
 L = @(t) interp1([0 0.2 0.3 0.4 tfin], ...
 [0 0 3e2 0 0 ], t, 'linear');
 M = @(t) 0;
 N = @(t) 0;
 INIT = [0; 0; 10; 0; 0; 0];
 case 4
 A = 1; B = 1; C = 10;
 L = @(t) 500*t;
 M = @(t) 0;
 N = @(t) 0;
 INIT = [0; 0; 50; 0; 0; 0];
end

Questa prima parte del codice è molto semplice. Si limita a pulire la memoria, definire due variabili del problema e mostrare un menù di scelta che presenta alcuni casi di problemi preimpostati. L’intero gruppo switch può essere ovviamente rimpiazzato da uno solo dei casi.
Notare che le leggi temporali delle componenti del momento L, M ed N, sono state definite utilizzando le anonymous function, che consistono nell’unico caso in cui è possibile definire una function all’interno di uno script in MATLAB. Queste funzioni sono state definite utilizzando un’altra comodissima routine, che è interp1: questa funzione prende in input (nel nostro caso) due vettori (uno di tempi e uno dei valori del momento) ed uno scalare. Utilizzando una semplice interpolazione lineare permette di restituire un valore della funzione al tempo t, interpolando i valori di M che sono stati forniti a determinati istanti di tempo. In pratica, tu gli dici quanto vale M al tempo 0, 1 e 5, lui ti dice quanto vale M in un qualsiasi istante compreso nell’intervallo [1,5], interpolando linearmente. Le funzionalità di interp1 sono molto più vaste, ma nel nostro caso la utilizziamo semplicemente per questo.

Vediamo ora come impostare il calcolo vero e proprio.

Mostra codice ▼

% Integrazione
opt = odeset('RelTol', prec, 'AbsTol', prec*ones(1,6));
hW = waitbar(0, 'Integrazione in corso...');
Par = {L M N A B C tfin hW};
[T, Y] = ode45(@eulero, [0 tfin], INIT, opt, Par);
delete(hW)

psi = Y(:,4);
theta = Y(:,5);
phi = Y(:,6);

La funzione odeset permette di impostare i parametri relativi all’algoritmo di integrazione. Senza addentrarci troppo nel suo funzionamento, ci limitiamo a settare i parametri di precisione, il cui valore influenza pesantemente il tempo di integrazione.

La funzione waitbar, come dice il nome, permette di creare una barra di caricamento a partire da un valore compreso tra 0 e 1, che nel nostro caso sarà l’intervallo di tempo di integrazione normalizzato. È una funzione comodissima, da quando l’ho scoperta la piazzo praticamente ovunque 😀

La funzione ode45 è quella che effettua l’integrazione. Per il funzionamento nel dettaglio si rimanda all’help di MATLAB, ma per quanto riguarda il nostro problema è sufficiente sapere che accetta in input una funzione (nel nostro caso eulero.m) che restituisce un vettore di sei righe, il cui valore è quello del secondo membro del sistema di equazioni differenziali (opportunamente modificato in modo da avere a primo membro le sole derivate prime delle variabili).

La funzione che andiamo ad integrare dipende dal tempo t, dal vettore di stato y e da una serie di parametri aggiuntivi Par:

Mostra codice ▼

function M = eulero(t, y, Par)
[L M N A B C tfin hW] = Par{:};

M = [(L(t) + y(2)*y(3)*(B-C))/A;
 (M(t) + y(1)*y(3)*(C-A))/B;
 (N(t) + y(1)*y(2)*(A-B))/C;
 y(1);
 y(2);
 y(3)];
 waitbar(t/tfin, hW)

Questa funzione viene richiamata ad ogni step di integrazione. Il suo compito è valutare le derivate prime del vettore di stato y. Oltre a fare ciò, questa funzione aggiorna anche la barra di caricamento ad ogni step. In questo modo potremo sapere graficamente a che punto dell’integrazione siamo arrivati.

La terza parte del codice, la più divertente, è quella che mostra graficamente i risultati dell’integrazione.

Mostra codice ▼

%% Grafica
% Diagramma angoli di Eulero
scrsz = get(0,'ScreenSize');
figure('position', [0 0 scrsz(3)/2 scrsz(4)])

subplot(3,1,[1 2])
plot(T, psi*deg2rad, ...
 T, mod(theta*deg2rad, 360) - 180, ...
 T, mod(phi*deg2rad, 360))
legend('\psi','\theta','\phi')
title('Angoli di Eulero')

subplot(3,1,3)
plot(T, L(T),'r-', T, M(T),'b--',T, N(T),'k-.')
legend('M_x','M_y','M_z')
title('Legge forzante')

% Animazione del moto
fH = figure('position', [scrsz(3)/2 0 scrsz(3)/2 scrsz(4)]);
% Mov = getframe;
Track = zeros(length(T), 3);

% Linearizzo il tempo
NtL = 300*tfin;
tL = linspace(0, tfin, NtL);
psi_L = interp1(T, psi, tL);
theta_L = interp1(T, theta, tL);
phi_L = interp1(T, phi, tL);

% Animazione
for i = 1:NtL
 figure(fH)
 clf, hold on, view(3), grid, camproj('perspective')
 title(sprintf('Tempo %.3g di %g', tL(i), tfin))
 % Matrice dei coseni direttori xyz
 dcm = ang2dcm(psi_L(i), theta_L(i), phi_L(i));
 % Asse
 Vdown = dcm * [0; 0; -4]; Vup = dcm * [0; 0; 4];
 plot3([Vup(1) Vdown(1)], [Vup(2) Vdown(2)], [Vup(3) Vdown(3)], ...
 'r', 'linewidth', 2)
 % Cilindro
 [Xcy, Ycy, Zcy] = cylinder(2, 30);
 CylD = dcm * [Xcy(1,:); Ycy(1,:); Zcy(1,:)];
 CylU = dcm * [Xcy(2,:); Ycy(2,:); Zcy(2,:)];
 surf([CylD(1,:); CylU(1,:)], ...
 [CylD(2,:); CylU(2,:)], ...
 [CylD(3,:); CylU(3,:)]);
 % Percorso dell'asse
 Track(i, 🙂 = Vup(:);
 if i <= 150, st = 1; else st = i-150; end
 plot3(Track(st:i,1), Track(st:i,2), Track(st:i,3), 'b')
 axis([-4 4 -4 4 -4 4])
 xlabel('x'), ylabel('y'), zlabel('z')
 drawnow;
 % Mov(i) = getframe(fH);
end

plot3(Track(:,1), Track(:,2), Track(:,3), 'b')

% movie2avi(Mov,'vudei.avi','fps',30)
% movie(Mov,30)

Vengono create due finestre: nella prima viene mostrato l’andamento degli angoli di Eulero, nella seconda l’animazione di un corpo che si muove secondo gli angoli di Eulero trovati.

La prima figura è molto semplice. L’unica particolarità è che gli angoli e vengono limitati agli intervalli [-180, 180] e [0, 360].

Prima di procedere con l’animazione è necessario conoscere un comportamento peculiare di ode45: insieme al vettore Y contenente le variabili di stato integrate, viene restituito anche un vettore temporale T cui corrispondono gli istanti di tempo relativi ad Y. Questo vettore è, generalmente, non lineare. ode45 ha infatti la capacità di modificare il passo di integrazione rendendolo più o meno fitto a seconda della variazione del vettore di stato. Se utilizzassimo questo vettore senza modificarlo, il risultato sarebbe un video che va piano quando ci sono variazioni irregolari di Y e veloce quando Y è più regolare.
Per ottenere un’animazione fisicamente realistica dobbiamo quindi linearizzare il tempo.

In questo caso, ho deciso di prendere 300 time-step per ogni secondo di simulazione (dovreste metterlo più grande se il vostro computer è performante, più piccolo nel caso contrario). La funzione linspace genera un vettore equidistanziato tra 0 e tfin. La funzione interp1 permette di valutare il vettore di stato in questi nuovi istanti di tempo, interpolando il vettore risultato di ode45.

Gli oggetti che vengono rappresentati graficamente sono tre: un cilindro, l’asse del disco e la sua traccia. Il cilindro ha sempre le stesse dimensioni, indipendentemente dai valori dei momenti di inerzia impostati. Ciò può essere fuorviante nel caso in cui A e B siano diversi tra loro (es. case 3), ma non importa.

Sia per la rappresentazione del cilindro che per l’asse, è necessario utilizzare la matrice dei coseni direttori. Questa matrice consente di rappresentare un generico vettore utilizzando gli angoli di Eulero. Per chi avesse installato l’Aerospace Toolbox di MATLAB, esiste una funzione angle2dcm che consente di scrivere una matrice di rotazione per una qualsiasi sequenza di rotazioni.
Dato che non tutti hanno questo strumento, ho aggiunto al codice la funzione ang2dcm, la quale non è altro che una versione semplificata di angle2dcm limitata alla sola sequenza XYZ, che è quella che ho utilizzato.

La traccia dell’asse viene mostrata soltanto per gli ultimi 150 timestep, in modo da non appesantire la rappresentazione grafica.

Nel codice sono presenti una serie di istruzioni commentate riferenti la variabile Mov e la funzione getframe. Questa funzione permette di acquisire l’immagine della figure per poterne creare successivamente un video. Se la utilizzate, tenete presente che richiede molta RAM e rallenta notevolmente il codice.

Assi giunti

Spesso i giroscopi si studiano utilizzando un sistema di riferimento diverso da quello utilizzato da me. Si considera il giroscopio montato su uno snodo cardanico (vedi figura) e si scrivono angoli e momenti rispetto agli assi degli snodi.

Questa formulazione è molto più conveniente da un punto di vista pratico, ma porta ad una formulazione delle equazioni più articolata. Senza addentrarci troppo nel discorso, la differenza sostanziale nell’utilizzo di ode45 per questo tipo di sistema è che una delle equazioni presenta due termini del secondo ordine. Per questo motivo è necessario impostare mediante odeset una matrice di massa (funzione come al solito di y e t) previa scrittura del sistema in questa forma:

La matrice di massa funziona esattamente come la funzione eulero.m. Notare che nel nostro caso la matrice di massa non è stata impostata, anche se questa corrisponderebbe alla matrice identità.

Video

Ecco infine i video dei 4 case proposti:

– Download del codice – eulero.zip

Il software xFoil è un potentissimo strumento che viene utilizzato in campo aeronautico per la progettazione dei profili alari subsonici. MATLAB è uno dei software più utilizzati nel campo dell’ingegneria. In questo articolo spiegherò come fare ad interfacciare i due software utilizzando la modalità di lavoro batch di xFoil.

Come appena detto, xFoil dà la possibilità di essere eseguito in modalità batch: un semplice file di testo contenente tutti i comandi in colonna può essere utilizzato come input del programma, che eseguirà il tutto in sequenza.
Il collegamento con MATLAB consiste nel creare il file di testo contenente i comandi, esportare i risultati di xFoil in dei file di testo e infine importare questi file all’interno di MATLAB.

I comandi principali che utilizzeremo per fare ciò sono:

• fopen e fprintf per creare un file di testo
• !xfoil < batch.dat per eseguire xFoil
• importdata per importare i risultati

Il codice

Scaricabile da questo link: Download script MATLAB-xfoil.zip
Mostra codice ▼

% Generazione e importazione di dati con il software xFoil
% Script elaborato da Alessandro Masullo, 03/12/2012
% Web: http://www.wirgilio.it
%
% Questo script crea un file di comandi per xFoil, lo esegue e importa i
% risultati nel workspace.
% Una cartella "dati" viene creata per includere i risultati di xFoil
% Un file di testo "comandi.dat" contiene i comandi che vengono fatti
% eseguire ad xFoil. In particolare viene esportato il diagramma del
% coefficiente di pressione di un NACA 4412 tra alfa -10 e alfa 15 con
% passo 0.5 a Reynolds 10^6
clc; clear;
%% Scrittura dei dati
% Crea una cartella per contenere i risultati di xFoil
if not(isdir('./dati')), mkdir('./dati'); end

% Comandi di xFoil
alfa = -10:0.5:15;
N = length(alfa);

comandi = {
 'NACA 4412'
 'OPER'
 'RE 1e6'
 'VISC ON'
 sprintf(['ALFA %.3f\n' ...
 'CPWR ./dati/cp(%d).dat\n'], [alfa; 1:N])
 'QUIT'
};

% Scrive i comandi nel file
id = fopen('comandi.dat', 'w+');
fprintf(id, '%s\n', comandi{:});
fclose(id);

% Esegue xfoil
!xfoil.exe < comandi.dat

%% Legge i risultati
% Apre in sequenza tutti i file generati da xfoil e li mette in un cell
% array
dati = cell(1, N);
for i = 1:N
 % La variabile foo è di appoggio e serve a importare il solo contenuto
 % numerico dei file generati da xfoil
 foo = importdata(sprintf('./dati/cp(%d).dat', i), ' ', 1);
 dati{i} = foo.data(:, : );
end

%% Post-processing
figure(1)
for i = 1:N
 plot(dati{i}(:, 1), dati{i}(:, 2), '-k')
 % Asse negativo verso il basso
 set(gca, 'ydir', 'reverse')
 % Calcola il cl da mettere nel titolo
 cl = trapz(dati{i}(:, 1), dati{i}(:, 2));
 title(sprintf('Diagramma del cp a cl = %.3f', cl))
 axis([0 1 -5 1])
 pause(0.5)
end

La spiegazione del codice

Per prima cosa, creiamo una variabile (cell array) contenente tutti i comandi da eseguire all’interno di xFoil. La variabile comandi è molto intuitiva e potrebbe essere utilizzata inserendo tutti i singoli comandi uno ad uno, ma la comodità di utilizzare MATLAB sta proprio nel non dover inserire manualmente i comandi ripetitivi, come ad esempio una sequenza di:

alfa 0
cpwr alfa(0).dat
alfa 0.1
cpwr alfa(1).dat
... e così via

Per cui, inseriamo manualmente i comandi principali, come la scelta del profilo e i parametri di viscosità, e utilizziamo la funzione sprintf per scrivere una sequenza di angoli di attacco.

comandi = {
 'NACA 4412'
 'OPER'
 'RE 1e6'
 'VISC ON'
};

La funzione sprintf, originaria del C, crea una stringa a partire da un modello contenente delle variabili. La sintassi è:

sprintf('stringa modello', var1, var2, ..., varN')

La struttura della stringa modello contiene una serie di “segnaposto”, identificati dal simbolo % (percento), all’interno del quale vengono inserite le variabili passate alla funzione come argomento. Es.:

sprintf('Il mio nome è %s, ho %d anni e sono alto %.2f m', 'Ale', 8, 1.753)
Il mio nome è Ale, ho 8 anni e sono alto 1.75 m

Il segnaposto %s sta per stringa, %d sta per numero intero, %.2f sta per numero decimale con 2 cifre dopo la virgola.

Nel nostro caso la stringa che vogliamo generare è:

alfa -10
cpwr cp(1).dat
alfa -9.5
cpwr cp(2).dat
...

Abbiamo quindi bisogno di un vettore contenente gli angoli di attacco (da -10 a 15 con passo 0.5)

alfa = -10:0.5:15 

E un vettore di indice (da 1… a quanti sono gli alfa)

N = length(alfa)
indice = 1:N 

Per praticità, utilizzeremo una caratteristica del modo di funzionare di sprintf su MATLAB. Se sprintf riceve come argomento una matrice, i valori che metterà nei segnaposto saranno gli elementi della matrice, presi in ordine procedendo lungo le colonne.
Dato che la nostra stringa vede in sequenza l’angolo d’attacco e il suo indice, potremo utilizzare una matrice di questo tipo:

[alfa(1), alfa(2), ..., alfa(n);
1, 2, ..., n]
Sarà interpretato come:
alfa(1), 1, alfa(2), 2, ... alfa(n), n

Per generare la sequenza di angoli d’attacco sarà quindi sufficiente utilizzare il comando:

 sprintf('ALFA %.3f\nCPWR ./dati/cp(%d).dat\n', [alfa; 1:N])

Dove il simbolo \n sta per “a capo”.

Da cui l’espressione completa dei comandi:

alfa = -10:0.5:15;
N = length(alfa);
comandi = {
 'NACA 4412'
 'OPER'
 'RE 1e6'
 'VISC ON'
 sprintf(['ALFA %.3f\n' ...
 'CPWR ./dati/cp(%d).dat\n'], [alfa; 1:N])
 'QUIT'
};

Questi comandi vanno inseriti in un file di testo, per cui creiamo un nuovo file e inseriamovi il contenuto della variabile comandi:

id = fopen('comandi.dat', 'w+');
fprintf(id, '%s\n', comandi{:});
fclose(id);

Infine, eseguiamo xFoil con questo comando:

!xfoil.exe < comandi.dat

Attenzione! Affiché questo comando venga eseguito correttamente, il file eseguibile xfoil.exe deve essere collocato nella vostra cartella di lavoro. Per conoscere la vostra cartella di lavoro digitate pwd nella command window.

Se è andato tutto a buon fine, nella cartella di lavoro dovreste trovarvi una nuova cartella, dati, contenente un bel po’ di file cp(x).dat. Sono ovviamente i file generati da xFoil.

A questo punto non ci resta che importare i risultati. Per fare ciò utilizzeremo la funzione importdata, la cui sintassi è:

importdata('nome del file', 'delimitatore', righe di intestazione)

Nel nostro caso il nome del file è una stringa del tipo cp(%d).dat, dove %d è un numero che va da 1 a quanti sono gli alfa.
Per importare tutti i file sarà sufficiente un ciclo for che, da 1 ad N, importa tutti i risultati di xFoil e li carica in un cell array di nome dati:

dati = cell(1, N);
for i = 1:N
 % La variabile foo è di appoggio e serve a importare il solo contenuto
 % numerico dei file generati da xfoil
 foo = importdata(sprintf('./dati/cp(%d).dat', i), ' ', 1);
 dati{i} = foo.data(:, : );
end

Come esplicitato nel commento, importdata cercherà di importare anche la prima riga di testo di ogni file, la quale ci ricorda che la prima colonna contiene le ascisse e la seconda i valori del Cp. Per evitare di importare tra i dati anche questa stringa di testo inseriamo il solo valore data all’interno dell’array di dati.

L’ultima parte del codice si limita a fare una serie di diagrammi in sequenza dei vari Cp calcolati da xFoil.

È possibile scaricare un archivio contenente questo script funzionante e la versione eseguibile di xFoil per Windows al seguente link:

Download script MATLAB-xfoil.zip

Variazione del numero di Reynolds

Ho ricevuto un’e-mail in cui mi si chiedeva una modifica del codice che permettesse la variazione del numero di Reynolds. Ho pensato che potesse essere utile anche a qualcun altro, e così l’ho pubblicata qui.
Non sarà certo il modo più elegante di farlo, ma le modifiche al codice sono minimizzate! La matrice matp, che viene letta da fprintf lungo le colonne, ha questa struttura:

Re(1)   Re(1)   Re(1)
alfa(1) alfa(2) alfa(3)
1       4       7
Re(2)   Re(2)   Re(2)
alfa(1) alfa(2) alfa(3)
2       5       8
Re(3)   Re(3)   Re(3)
alfa(1) alfa(2) alfa(3)
3       6       9

Mostra codice ▼

clc; clear;
%% Scrittura dei dati
% Crea una cartella per contenere i risultati di xFoil
if not(isdir('./dati')), mkdir('./dati'); end

% Comandi di xFoil
alfa = 0:15;
Re = [1e5 5e5 1e6 5e6 1e7 5e7 1e8 1e9];
Na = length(alfa);
NRe = length(Re);
N = NRe*Na;

matp = zeros(NRe*3, Na);
matp(1:3:end-2, : ) = repmat(Re', 1, Na);
matp(2:3:end-1, : ) = repmat(alfa, NRe, 1);
matp(3:3:end, : ) = reshape(1:N, NRe, Na);

comandi = {
 'NACA 4412'
 'OPER'
 'ITER 200'
 'VISC'
 '1e6'
 sprintf(['RE %g\n' ...
 'ALFA %.3f\n' ...
 'CPWR ./dati/cp(%d).dat\n' ...
 'INIT\n'], matp)
 'QUIT'
};

% Scrive i comandi nel file
id = fopen('comandi.dat', 'w+');
fprintf(id, '%s\n', comandi{:});
fclose(id);

% Esegue xfoil
!xfoil < comandi.dat

%% Legge i risultati
% Apre in sequenza tutti i file generati da xfoil e li mette in un cell
% array
dati = cell(1, N);
for i = 1:N
 % La variabile matp è di appoggio e serve a importare il solo contenuto
 % numerico dei file generati da xfoil
 foo = importdata(sprintf('./dati/cp(%d).dat', i), ' ', 1);
 dati{i} = foo.data(:, : );
end

%% Post-processing
figure(1)
for i = 1:N
 plot(dati{i}(:, 1), dati{i}(:, 2), '-k')
 % Asse negativo verso il basso
 set(gca, 'ydir', 'reverse')
 title(sprintf('Diagramma del cp ad alfa = %.1f, Re = %g,', matp(3*i-1), matp(3*i-2)))
 axis([0 1 -5 1])
 pause(0.5)
end

Suggerimenti

Ai fini delle prestazioni, xFoil risulta estremamente più veloce se sopprimete l’output grafico, che vi sarà probabilmente inutile se avete intenzione di lavorare con MATLAB sui risultati.
Per fare ciò, potete inserire questi tre elementi ai primi comandi della variabile comandi:

comandi = {
'plop'
'g'
''
% tutto il resto

}

Accertatevi inoltre, specialmente se disabilitate l’output grafico, che nei casi estremi che andate a studiare xFoil giunga a convergenza, altrimenti importerete dei risultati che potrebbero essere sballati!
Scegliete quindi un numero di iterazioni opportunamente elevato per assicurarvi di ciò.

Troppo spesso mi è capitato di imparare ad utilizzare MATLAB e poi di dimenticare tutto a causa di lunghi periodi di inutilizzo.

“Come si chiamava quella funzione?”
“Com’è che si faceva a fare… ?”
“Che sintassi usavo per fare quel coso… ?”

Ma stavolta ho detto basta! E così mi sono deciso a scrivere un file in cui ho riassunto tutte le funzioni, i costrutti e le cose (imho) più utili da ricordare di MATLAB. Il risultato è questo pratico cheatsheet  che ho deciso di condividere con tutti.

Se avete consigli o suggerimenti per migliorarlo, scrivetelo nei commenti.

Download MATLAB cheatsheet (.PDF) – 92 kB