Transcript Studio di sistemi di equazioni differenziali lineari del primo ordine

Modelli Dinamici
(anche per le questioni di coppia)
Studio di sistemi di equazioni differenziali lineari
del primo ordine
R
Catena Tiziana
J
Strogatz nel 1994 affronta con una divertente metafora lo
studio di sistemi di equazioni differenziali ordinarie del primo
ordine in forma esplicita. Immagina che le due funzioni, R(t) e
J(t) rappresentino l'amore che Romeo e Giulietta provano
l'uno per l'altro, in ogni istante t. Il nostro obiettivo è di capire
come variano i sentimenti dei due protagonisti a seconda delle
condizioni iniziali e dei parametri a,b,c,d.
a e d rappresentano quanto il soggetto sia motivato dall'amore
che prova per l'altro, c e d rappresentano invece la sensibilità
verso il sentimento provato dall'altro per loro. Tutti i parametri,
come le condizioni iniziali, possono avere segno positivo e
negativo
Il sistema in considerazione è autonomo, cioè le sue
equazioni non dipendono esplicitamente dalla variabile
indipendente t. In ogni momento ogni funzione dipende anche
dal valore che l'altra sta assumendo in quell'istante.
Stabilità del sistema differenziale
Un sistema di equazioni differenziali è stabile se TUTTI gli
autovalori dello Jacobiano hanno parte reale <0. Vediamo perché.
Il nostro sistema può essere rappresentato come
J
Notiamo che J corrisponde proprio alla matrice delle derivate
parziali delle due derivate (cioè lo Jacobiano).
Se riusciamo a diagonalizzare la matrice J possiamo sostituire il
sistema con uno equivalente più semplice.
J
Che si traduce
nel sistema
Ma questo sistema ha soluzione:
Se quindi la parte reale degli autovalori è
positiva, al crescere della variabile t le due
funzioni tendono a valori infiniti il cui valore
è impossibile da determinare perché
instabile. Se invece abbiamo parte reale
negativa, le funzioni tenderanno ad un certo
valore che si mantiene
approssimativamente costante per
.
Esempio di sistema non stabile:
heun_sistema(-2,-20,4,5)
Autovalori: 3.5 + 8.8176i
3.5 – 8.8176i
Condizioni sui coefficienti
Abbiamo dunque visto che la parte reale degli
autovalori dello Jacobiano del sistema devono
essere minori di zero. Calcoliamo gli autovalori di J
in funzione dei parametri e vediamo come le
condizioni di stabilità si ripercuotono sui coefficienti
a,b,c,d.
1,2
A questo punto abbiamo due casi,
Discutiamoli.
oppure
.
: in tal caso gli autovalori saranno complessi e la parte reale
sarà (a+d)/2. Il sistema allora sarà ben condizionato se:
: in questo caso dobbiamo assicurarci che ENTRAMBI gli
autovalori siano negativi. Per l'autovalore in cui si considera il
segno + davanti alla radice dobbiamo assicurarci non solo che
(a+d)<0 ma che il suo modulo sia anche maggiore di quello della
radice. Invece, per l'autovalore che ha segno – davanti alla radice
è accettabile sia che (a+d)<0 sia che (a+d)>0, assicurandoci però
che il suo modulo sia minore del contributo negativo della radice.
Ma dato che le condizioni di ENTRAMBI gli autovalori devono
essere soddisfatte abbiamo come risultato totale le seguenti
limitazioni:
Metodi a confronto
Per l'analisi grafica del sistema differenziale utilizzeremo tre
metodi iterativi;
Metodo di Runge-Kutta classico
Metodo di Heun
Metodo di Eulero
Tutti e tre i metodi sono consistenti e stabili, quindi convergenti.
Tutti e tre approssimano la curva calcolando prima la tangente
alla curva nel punto iniziale e poi costruendo successive
approssimazioni sulla tangente costruita al passo precedente.
Tuttavia il metodo di Runge-Kutta che sfrutta la previsione delle
derivate successive, risulta essere il più preciso, seguito dal
metodo di Heun. Vediamo qualche esempio pratico che mette
a confronto l'efficienza dei diversi metodi e l'influenza della
scelta di diversi passi di integrazione h.
Nota: I codici delle rispettive function elaborate su Octave si trovano a fine
presentazione.
Metodi a confronto
Eulero, h=0.5
Eulero, h=0.05
Runge-Kutta, h=0.5
Runge-Kutta, h=0.05
Metodo di Heun con diversi h
h=0.03
h=0.04
Anche se inizialmente le due
curve sono molto simili,
andando avanti con le
iterazioni l'errore insito nella
curva rossa si amplifica
provocando perdita di
informazioni.
Cambiando i parametri del
sistema...
Rimanendo sulla metafora proposta da Strogatz, proviamo ora a
cambiare i valori di a,b,c,d per vedere come evolve l'amore dei
due protagonisti.
In un grafico rappresenteremo l'andamento delle due funzioni nel
tempo, nell'altro si rappresenterà lo spazio delle fasi, ovvero lo
stato di un'equazione in funzione dell'altra.
Romeo sarà rappresentato dalla curva rossa e Giulietta dalla
curva blu.
Per il primo grafico facciamo arrestare la funzione ad un certo
tempo t da noi stabilito, per il secondo grafico la facciamo operare
finché non trova i valori di equilibrio, così da avere l'andamento
complessivo del sistema differenziale.
Controlliamo sempre se i parametri rispettano le condizioni di
stabilità e se effettivamente gli autovalori dello Jacobiano hanno
parte reale <0.
La coppia facilmente condizionabile
Supponiamo di porre a=d=0, e quindi di immaginare che R e J
non vengano affatto motivati dal proprio amore per l'altro, ma
che tutto sia determinato dal comportamento del compagno.
Partendo da una condizione iniziale positiva con R(1)=1 e J(1)=1
avremo:
Con b>0,c>0: situazione di amore eterno sempre crescente.
(Diventa invece un odio crescente se le condizioni iniziali hanno
valori negativi)
Con b<0,c<0: uno ama eternamente, l'altro odia terribilmente.
Con b e c discordi: ciclo infinito di amore ed odio che non
termina mai, più il tempo passa più sia l'odio che l'amore sono
forti.
a=d=0 → b>0,c>0 e R(1),J(1) positivi
a=d=0 → b>0,c>0 e R(1),J(1) negativi
a=d=0 → b<0,c<0
a=d=0 → b*c<0
Se Romeo 'fa il duro'
Supponiamo che l'amore di Romeo sia completamente insensibile ai
cambiamenti sentimentali della consorte. La funzione R sarà quindi
una costante, con a=0 e b=0. In questo caso l'evoluzione
dell'amore dipenderà sia dal segno dell'amore che Romeo prova
costantemente, sia dal carattere di Giulietta.
Se c>0,d>0 (J crede nel suo amore e in quello del compagno):
Se Romeo prova qualcosa per lei, il suo amore crescerà sempre di
più, se invece lui non la ama lei lo odierà per tutta la vita.
Se c>0,d<0 (J crede nell'amore di R ma non nel suo):
in questo caso R renderà costanti anche i sentimenti della consorte.
Se R prova qualcosa per lei, l'amore di J scenderà fino ad un certo
valore comunque positivo, poi rimarrà costante per sempre, se
invece R non la ama, J comincerà ad odiarlo ma oltre un certo valore
il suo odio diventerà costante. Comunque, non ci sono oscillazioni.
a=b=0 → c>0,d>0
R(1)>0
R(1)>0
R(1)<0
R(1)<0
R(1)>0
a=b=0 → c>0,d<0
R(1)>0
R(1)<0
R(1)<0
La coppia degli opposti
È vero che gli opposti si attraggono? Beh, a sentire la
matematica no. Supponiamo infatti di prendere dei
parametri tali che J e R reagiscano in maniera
opposta ai propri sentimenti e ai sentimenti dell'altro.
Avremo c e b, e d e a con segni opposti.
Quello che si verifica, a prescindere dalle condizioni
iniziali, è un iniziale alternarsi caotico di amore e odio
che tende ad affievolirsi con il tempo, fino a
raggiungere la completa apatia. (È l'opposto di
quanto accadeva sempre per b*c<0 ma con a=d=0).
c,b e d,a hanno segno opposto
Stati stazionari
Ci proponiamo infine di analizzare quand'è che
questo sistema differenziale raggiunge la
stazionarietà, ovvero quando le soluzioni rimangono
costanti. Graficamente abbiamo già notato qualcosa,
arriviamo però a una condizione più generale per via
analitica.
Trovare i valori stazionari del sistema significa
trovare il vettore soluzione in cui entrambe le derivate
prime si annullano. Ma dato che conosciamo la forma
delle derivate prime trovare i punti in cui si annullano
corrisponde a risolvere un sistema lineare del tipo:
Dobbiamo quindi trovare il nucleo della trasformazione che ha
per matrice la matrice dei coefficienti.
Abbiamo considerato però che un sistema di equazioni
differenziali è ben condizionato solamente quando la parte
reale di TUTTI gli autovalori dello Jacobiano è <0, allora
stiamo supponendo che non esistano autovalori nulli
nemmeno per questa matrice che corrisponde proprio allo
Jacobiano del sistema. Se non esistono autovalori nulli allora
la matrice è sicuramente invertibile, e di conseguenza l'unico
vettore appartenente al nucleo è il vettore nullo. Non a caso
infatti, molti degli esempi precedenti vedevano le funzioni
convergere allo zero appiattendosi. Eccone un altro esempio:
Se però consideriamo anche sistemi non perfettamente
ben condizionati, ma cerchiamo punti stazionari anche in
quei sistemi che si avvicinano molto alle condizioni di
stabilità, allora il vettore nullo, cioè il punto (0,0), non è
l'unico punto stazionario. Troviamo ad esempio le
soluzioni del SLO:
Lo spazio delle soluzioni comprende tutti i vettori soluzione
di forma
per i quali le derivate prime del sistema si
annullano. Ogni vettore soluzione ha coordinate multiple di
quelle dell'autovettore base.
Non possiamo però fare infiniti tentativi
per trovare il valore di t o di k che
effettivamente ci restituisce valori di R e J
appartenenti allo spazio delle soluzioni e
quindi multipli delle componenti del
vettore base. Andiamo perciò a
controllare che il rapporto tra le
coordinate di R e J nei punti in cui
sembrano appiattirsi, è uguale al rapporto
tra le componenti dell'autovettore.
I punti stazionari del sistema saranno
proprio quei punti in cui il rapporto tra
R e J coincide con quello previsto.
Piccola nota: a conferma dell'instabilità del sistema si è
dovuta interrompere la funzione a causa di un crash di
Octave.
>>> heun_sistema(2,-10,4,-20)
Function “eulero_sistema”
function[]=eulero_sistema(a,b,c,d)
%Metodo di Eulero per la soluzione di sistemi a due dimensioni
%di equazioni differenziali di primo ordine omogenee.
%Inserire i coefficienti sistema di equazioni differenziali.
printf('dy= @(x,y) a.*x+b.*y\ndx= @(x,y) c.*x+d.*y\n\n');
dx= @(x,y) a.*x+b.*y;
dy= @(x,y) c.*x+d.*y;
x(1)=input('Inserire valore → x(1)=');
y(1)=input('Inserire valore → y(1)=');
t(1)=1;
h=input('Inserire il passo di integrazione → h=');
if isempty(h) h=0.1;
printf('\nAttenzione, valore non inserito!\n\n') endif
r=input('Per quale t si vuole calcolare la soluzione del sistema?:');
if isempty(r) r=10;
printf('\nAttenzione, valore non inserito!\n\n') endif
err=input('Precisione con cui calcolare valori di equilibrio del sistema → err=');
if isempty(err) err=0.5*10^-5; endif
%------------------------------------------ricerca della soluzione e plot fino a quel punto
printf('\n t
t(n)
y(n)
x(n)\n');
printf(' %d\t
%f\t
%f\t
%f\t\n', [1 t(1) y(1) x(1)]);
for (n=2:r/h)
%r/h è il numero delle iterazioni per raggiungere quel
tempo
t(n)=(t(1)+(n-1)*h);
%vettore tempo
X=x(n-1);
%salvo i valori di x e y precedenti
Y=y(n-1);
x(n)=x(n-1)+dx(X,Y)*h; %calcolo le successive componenti di x e y
y(n)=y(n-1)+dy(X,Y)*h;
printf('
%d\t
%f\t
%f\t
%f\t\n', [n t(n) y(n)
x(n)]);
endfor
figure(1)%-------------------subplot(1,2,1);
plot(t,x,'r--*');
del tempo
hold on
plot(t,y,'b--*');
del tempo
title('Sistema fino al valore richiesto');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
%disegno vettore x in funzione
%disegno vettore y in funzione
figure(2)%------------------subplot(1,2,1)
plot(x,y,'k--*');
title('Spazio delle fasi fino al valore richiesto');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%-----------------------------------------------------------------ricerca del punto di
equilibrio
n=1;
err_ass=10;
%inizializzazione casuale per
entrare nel ciclo
while(err_ass>err)
n=n+1;
t(n)=(t(1)+(n-1)*h);
X=x(n-1);
Y=y(n-1);
x(n)=x(n-1)+dx(X,Y)*h;
y(n)=y(n-1)+dy(X,Y)*h;
err_ass=abs(y(n)-x(n));
endwhile
printf(' --------\n');
printf('\n
↓ Valori di equilibrio del sistema (err= %f) ↓\n',[err]);
printf(' %d\t
%f\t
%f\t\n', [t(n) y(n) x(n)]);
figure(1)%--------------------subplot(1,2,2)
plot(t,x,'r--*');
hold on
plot(t,y,'b--*');
title('Sistema fino ai valori di equilibrio');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
figure(2)%---------------------subplot(1,2,2)
plot(x,y,'k--*');
title('Spazio delle fasi fino al valore di equilibrio');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%------------------------------------------------------------------------------------------------
Function “heun_sistema”
function[]=heun_sistema(a,b,c,d)
%Metodo di Heun per la soluzione di sistemi a due dimensioni
%di equazioni differenziali di primo ordine omogenee.
%Inserire i coefficienti sistema di equazioni differenziali.
%printf('dy= @(x,y) a.*x+b.*y\ndx= @(x,y) c.*x+d.*y\n\n');
dx= @(x,y) a.*x+b.*y;
dy= @(x,y) c.*x+d.*y;
x(1)=input('Inserire valore → x(1)=');
y(1)=input('Inserire valore → y(1)=');
t(1)=0;
h=input('Inserire il passo di integrazione → h=');
if isempty(h) h=0.1;
printf('\nAttenzione, valore non inserito!\n\n') endif
r=input('Per quale t si vuole calcolare la soluzione del sistema?:');
if isempty(r) r=10;
printf('\nAttenzione, valore non inserito!\n\n') endif
err=input('Precisione con cui calcolare valori di equilibrio del sistema → err=');
if isempty(err) err=0.5*10^-5; endif
%---------------------------------------ricerca della soluzione e plot fino a quel punto
printf('\iterazione n
t(n)
y(n)
x(n)\n');
printf('%d\t
%f\t
%f\t
%f\t\n', [1 t(1)
y(1) x(1)]);
for(n=2:(r/h)+1)
t(n)=(t(1)+(n-1)*h);
Y=y(n-1);
X=x(n-1);
%decido numero di iterazioni
%vettore tempo
y(n)=Y + (h/2)*( dy(X,Y) + dy(X+h*dx(X,Y), Y+h*dy(X,Y)) );
x(n)= X + (h/2)*( dx(X,Y) + dx(X+h*dx(X,Y), Y+h*dy(X,Y)) );
printf('%d\t
%f\t
%f\t
%f\t\n', [n t(n) y(n) x(n)]);
endfor
figure(1)%-------------------subplot(1,2,1);
plot(t,x,'r--*');
funzione del tempo
hold on
plot(t,y,'b--*');
funzione del tempo
title('Sistema fino al valore richiesto');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
figure(2)%------------------subplot(1,2,1)
plot(x,y,'k--*');
title('Spazio delle fasi fino al valore richiesto');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%-----------------------------------------------------------------ricerca del punto di
equilibrio
n=1;
err_ass=10;
%inizializzazione casuale per
entrare nel ciclo
while(err_ass>err)
n=n+1;
t(n)=(t(1)+(n-1)*h);
X=x(n-1);
Y=y(n-1);
y(n)=Y + (h/2)*( dy(X,Y) + dy(X+h*dx(X,Y), Y+h*dy(X,Y)) );
x(n)= X + (h/2)*( dx(X,Y) + dx(X+h*dx(X,Y), Y+h*dy(X,Y)) );
err_ass=abs(y(n)-x(n));
endwhile
printf(' --------\n');
printf('\n
↓ Valori di equilibrio del sistema (err= %f) ↓\n',[err]);
printf('
%d\t
%f\t
%f\t\n', [t(n) y(n) x(n)]);
figure(1)%--------------------subplot(1,2,2)
plot(t,x,'r--*');
hold on
plot(t,y,'b--*');
title('Sistema fino ai valori di equilibrio');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
%disegno vettore x in
%disegno vettore y in
figure(2)%---------------------subplot(1,2,2)
plot(x,y,'k--*');
title('Spazio delle fasi fino al valore di equilibrio');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%--------------------------------------------------------------------------------------------------
Function “runge_kutta_sistema”
function[]= runge_kutta_sistema(a,b,c,d)
%Metodo di Runge_Kutta per la soluzione di sistemi a due dimensioni
%di equazioni differenziali di primo ordine omogenee.
%Inserire i coefficienti sistema di equazioni differenziali.
%printf('dy= @(x,y) a.*x+b.*y\ndx= @(x,y) c.*x+d.*y\n\n');
dx= @(x,y) a.*x+b.*y;
dy= @(x,y) c.*x+d.*y;
x(1)=input('Inserire valore → x(1)=');
y(1)=input('Inserire valore → y(1)=');
t(1)=0;
h=input('Inserire il passo di integrazione → h=');
if isempty(h) h=0.1;
printf('\nAttenzione, valore non inserito!\n\n') endif
r=input('Per quale t si vuole calcolare la soluzione del sistema?:');
if isempty(r) r=10;
printf('\nAttenzione, valore non inserito!\n\n') endif
err=input('Precisione con cui calcolare valori di equilibrio del sistema → err=');
if isempty(err) err=0.5*10^-5; endif
%--------------------------------------------------ricerca della soluzione e plot fino a quel punto
printf('\n iterazione n
t(n)
y(n)
x(n)\n');
printf('%d\t
%10.5f\t
%10.5f\t
%10.5f\t \n', [1 t(1) y(1) x(1)]);
for(n=2:(r/h))
%decido numero di iterazioni
t(n)=(t(1)+(n-1)*h);
Y=y(n-1);
X=x(n-1);
%vettore tempo, ogni componente è incrementata di h
K1= @(x,y) dx ( x , y );
k1= @(x,y) dy ( x , y );
K2=@(x,y) dx ( x+(h/2)*K1(x,y) , y+(h/2)*k1(x,y) );
k2=@(x,y) dy ( x+(h/2)*K1(x,y) , y+(h/2)*k1(x,y) );
K3=@(x,y) dx ( x+(h/2)*K2(x,y) , y+(h/2)*k2(x,y) );
k3=@(x,y) dy ( x+(h/2)*K2(x,y) , y+(h/2)*k2(x,y) );
K4=@(x,y) dx ( x+h*K3(x,y) , y+h*k3(x,y) );
k4=@(x,y) dy ( x+h*K3(x,y) , y+h*k3(x,y) );
y(n)= Y+(h/6)*( k1(X,Y) + 2*k2(X,Y) + 2*k3(X,Y) + k4(X,Y));
x(n)= X+(h/6)*( K1(X,Y) + 2*K2(X,Y) + 2*K3(X,Y) + K4(X,Y));
printf('%1.5d\t
%10.5f\t %10.5f\t
%10.5f\t\n', [n t(n) y(n) x(n)]);
endfor
figure(1)%-------------------subplot(1,2,1);
plot(t,x,'r--*');
%disegno vettore x in funzione del tempo
hold on
plot(t,y,'b--*');
%disegno vettore y in funzione del tempo
title('Sistema fino al tempo richiesto');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
figure(2)%------------------subplot(1,2,1)
plot(x,y,'k--*');
title('Spazio delle fasi fino al tempo richiesto');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%----------------------------------------------------------------------------ricerca del punto di equilibrio
n=1;
err_ass=10;
%inizializzazione casuale per entrare nel ciclo
while(err_ass>err)
n=n+1;
t(n)=(t(1)+(n-1)*h);
%vettore tempo, ogni componente è incrementata di h
Y=y(n-1);
X=x(n-1);
K1= @(x,y) dx ( x , y );
k1= @(x,y) dy ( x , y );
K2=@(x,y) dx ( x+(h/2)*K1(x,y) , y+(h/2)*k1(x,y) );
k2=@(x,y) dy ( x+(h/2)*K1(x,y) , y+(h/2)*k1(x,y) );
K3=@(x,y) dx ( x+(h/2)*K2(x,y) , y+(h/2)*k2(x,y) );
k3=@(x,y) dy ( x+(h/2)*K2(x,y) , y+(h/2)*k2(x,y) );
K4=@(x,y) dx ( x+h*K3(x,y) , y+h*k3(x,y) );
k4=@(x,y) dy ( x+h*K3(x,y) , y+h*k3(x,y) );
y(n)= Y+(h/6)*( k1(X,Y) + 2*k2(X,Y) + 2*k3(X,Y) + k4(X,Y));
x(n)= X+(h/6)*( K1(X,Y) + 2*K2(X,Y) + 2*K3(X,Y) + K4(X,Y));
err_ass=abs(y(n)-x(n));
endwhile
printf('
---------\n');
printf('\n
↓ Valori di equilibrio del sistema (err= %f) ↓\n',[err]);
printf('%1.5d\t
%10.5f\t %10.5f\t
%10.5f\t\n', [n t(n) y(n)
x(n)]);
figure(1)%--------------------subplot(1,2,2)
plot(t,x,'r--*');
hold on
plot(t,y,'b--*');
title('Sistema fino ai valori di equilibrio');
xlabel('tempo(t)')
ylabel('x/y')
grid
hold off
figure(2)%---------------------subplot(1,2,2)
plot(x,y,'k--*');
hold on
plot(x(n),y(n),'*r')
title('Spazio delle fasi fino al valore di equilibrio');
xlabel('x(n)')
ylabel('y(n)')
grid
hold off
%--------------------------------------------------------------------------------------------------
E per finire…
Animazione dello spazio delle fasi con a=d=0, b*c<0
J
R
R
J