NOTA! Questo sito utilizza i cookie e tecnologie simili.

Se non si modificano le impostazioni del browser, l'utente accetta. Per saperne di piu'

Approvo

Matlab per la risoluzione delle equazioni differenziali (ODE), con esempio

 

Per chi ha la necessità di risolvere equazioni differenziali, ODE e loro sistemi, Matlab mette a disposizione alcune funzioni di libreria davvero utili, a seconda delle esigenze:

Ode45 Problemi non-stiff, ordine di accuratezza medio.
Ode23 Problemi non-stiff, ordine di accuratezza basso.
Ode113 Problemi non-stiff, ordine di accuratezza variabile (basso-alto).
Ode15s Problemi stiff, ordine di accuratezza variabile (basso-medio)
Ode23s Problemi stiff, ordine di accuratezza basso.
Ode23t Problemi "moderatamente stiff", ordine di accuratezza basso.
Ode23tb Problemi stiff, ordine di accuratezza basso.

 

Per utilizzare tali funzioni, è necessario definire l'equazione differenziale (o le equazioni differenziali, nel caso di un sistema) in un file .m, ove andranno definiti eventualmente anche altri parametri necessari per il calcolo, mentre bisognerà predisporre, nell'area di lavoro di Matlab, un array per l'output.

 

L'esempio seguente mostra come impostare il sistema di ODE per il problema di Keplero modificato (un tipo particolare di sistema hamiltoniano), composto da quattro equazioni con quattro incognite.

Ecco l'equazione del sistema hamiltoniano H:

ModifiedKeplerProblemHfunction

dove r è la radice quadrata della somma dei quadrati di q1 e q2. Il valore alpha va impostato a piacere; nel nostro caso, abbiamo scelto 0.01.

Le quattro incognite sono quindi p1, p2, q1, q2; scrivendo ModifiedKeplerProblemPvector ci riferiremo al "vettore" di variabili p, ossia in questo caso p1 e p2; discorso analogo va fatto per q.

 

Il sistema di ODE per tale problema è il seguente:

ModifiedKeplerProblemODEsystem

In questo caso, in un file .m a parte (funzioneKeplero.m, in questo caso) si crea l'array contenente le equazioni del sistema; il contenuto di tale file è il seguente:

%File funzioneKeplero.m
function funzione = funzioneKeplero(t, parametri) %parametri(1) = p1iniziale; parametri(3) = q1iniziale
alpha = 0.01; % Non posso specificarlo fuori e passarlo come quinto parametro: Matlab accetta solo un array di lunghezza 4 come secondo argomento per ode45
r = sqrt(parametri(3).^2 + parametri(4).^2);
coeff = - ((1 ./ (r.^3)) + ((3 .* alpha) ./ (2 .* r.^5)));
funzione = [coeff .* parametri(3); coeff .* parametri(4); parametri(1); parametri(2)]; %eq (nell'ordine) q1, q2, p1, p2

 

Adesso, nella Command Window di Matlab, bisogna fare in modo che una funzione della suite ODE del programma elabori, dato l'intervallo di tempo e il passo, tale sistema, depositando l'output (i valori delle quattro variabili p1, p2, q1, q2) nell'array "output", appunto; in realtà, output diverrà una matrice: ad ogni riga, si avranno i valori delle quattro variabili per un dato istante di tempo (ricordiamoci che, dato ad esempio l'intervallo [0, 500], impostando il passo h = 0.1 avremo, in realtà, (500-0) * (1/h) = 5000 righe-istanti di tempo, per cui richiamando la riga 3100 di output, per esempio, staremo esaminando l'istante temporale 310, e così via).

 

Ecco uno script che potete copiare nella Command Window di Matlab (ovviamente, dovrete avere il file funzioneKeplero.m nella vostra area di lavoro) per ottenere il grafico, come vedrete a video, il "diagramma di fase" delle variabili q1 e q2 per tale sistema nell'intervallo [0,500] con passo 0.1 con il metodo MidPoint (ODE23) e con Runge-Kutta del quarto ordine (ODE45).

>> %Istruzioni da eseguire nella Command Window di Matlab
>> a = 0;
>> b = 500;
>> h = 0.1;
>> intervallo = a:h:b;

>> beta = 0.6;
>> q1iniziale = 1- beta;
>> q2iniziale = 0;
>> p1iniziale = 0;
>> p2iniziale = sqrt((1+beta)/(1-beta));

>> parametri = [p1iniziale, p2iniziale, q1iniziale, q2iniziale];

>> %MidPoint (ODE23) con opzioni di tolleranza di default
>> [t, output] = ode23('funzioneKeplero', intervallo, parametri);

>> %RK4 (ODE45) con opzioni tolleranza impostate esplicitamente mediante ODESET
>> opzioniTolleranza = ODESET('RelTol', 1e-4, 'AbsTol', 1e-6); %Valori di default: 1e-3 e 1e-6, ossia 0.01% e 0.00001%
>> [t, output2] = ode45('funzioneKeplero', intervallo, parametri, opzioniTolleranza);

>> subplot(1,2,1);
>> xlabel('MidPoint esplicito');
>> plot(output(:,3), output(:,4));
>> subplot(1,2,2);
>> xlabel('RK4');
>> plot(output2(:,3), output2(:,4));
>>

 

Ecco l'immagine che otterrete a video:

OutputVideo-q1q2-Kepler

 

Non vi sarà sfuggita la riga di codice:

>> opzioniTolleranza = ODESET('RelTol', 1e-4, 'AbsTol', 1e-6); %Valori di default: 1e-3 e 1e-6, ossia 0.01% e 0.00001%

utilizzata per impostare i valori di tolleranza relativa ed assoluta per la funzione ODE da utilizzare.

 

In Matlab infatti è possibile definire due tipi di "tolleranze" nella risoluzione delle ODE: relativa ed assoluta. La relazione tra tali valori e l'errore ad ogni passo è la seguente (dalla documentazione di Matlab):

e(i) <= max(RelTol*abs(y(i)),AbsTol(i))

 

Per impostare valori differenti per questi parametri, bisogna quindi fare ricorso a ODESET.

 

ODESET consente di impostare DIVERSI parametri per la risoluzione delle ODE (per una trattazione più approfondita, si rimanda all'help di Matlab); nel nostro caso, abbiamo impostato solo i valori di tolleranza relativa ed assoluta da utilizzare con ODE45 (bisogna passare le impostazioni definite con odeset alle funzioni della libreria ODE di Matlab come quarto parametro).

 

 

COMMENTI, DOMANDE, ALTRE RISORSE, ...

Per commenti, domande o altro, accedere al Forum del sito (cliccare qui); per informazioni sulla Registrazione al sito, cliccare qui.

 

 

 
Vai all'inizio della pagina