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:

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:

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

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

Il sistema di ODE per tale problema è il seguente:

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 p₁, p₂, q₁, q₂) 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 q₁ e q₂ 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:

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).

*          *          *