Pràctica: DETERMINACIÓ D'ESDEVENIMENTS EN EDOs
En aquesta pràctica veurem com calcular resultats associats a una EDO o sistema d'EDOs, com poden ser el període d'una òrbita periòdica, la màxima amplitud d'una oscil·lació, etc. Això comporta resoldre equacions (és a dir, trobar zeros) en què la funció involucrada és la solució de l'EDO.
Una primera possiblitat seria usar un mètode genèric de càlcul de zeros, com és el mètode de Newton, aplicant-lo a la solució de l'EDO obtinguda amb ode45. Veurem però que podem obtenir el resultat de manera més directa utilitzant ode45 amb l'opció Events.
Exemple de sistema d'EDOs amb òrbites periòdiques
Com a il·lustració veurem com calcular el període d'una òrbita periòdica del sistema d'EDOs (autònom)
que té la funció 
com a quantitat conservada, i per tant les seves solucions recorren corbes de nivell d'aquesta funció. A més, el sistema té els punts d'equilibri 
i 
que són al mateix temps els punts crítics de la funció V (un mínim i un punt de sella). Representem el camp vectorial associat al sistema i, a continuació, les corbes de nivell de 
com a quantitat conservada, i per tant les seves solucions recorren corbes de nivell d'aquesta funció. A més, el sistema té els punts d'equilibri i que són al mateix temps els punts crítics de la funció V (un mínim i un punt de sella). Representem el camp vectorial associat al sistema i, a continuació, les corbes de nivell de close all %tanquem figures anteriors (optatiu)
figure %nova figura
hold on
P=@(x,y) y;
Q=@(x,y) -x+x.^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
x=-0.7:0.2:1.5; %finestra a visualitzar
y=-0.8:0.2:0.8;
[xg,yg]=meshgrid(x,y); %xarxa de (pocs) punts
Pg=P(xg,yg); %avaluem el camp vectorial sobre la xarxa
Qg=Q(xg,yg);
quiver(xg,yg,Pg,Qg,1.5,'Color','g') %camp vectorial per un factor
plot(0,0,'r*') %marquem els punts equilibri
plot(1,0,'r*')
axis equal
figure %nova figura
hold on
V=@(x,y) y.^2/2+x.^2/2-x.^3/3;
x=-0.7:0.01:1.5;
y=-0.8:0.01:0.8;
[xg,yg]=meshgrid(x,y); %nova xarxa de punts (mes fina)
zg=V(xg,yg); %avaluem la quantitat conservada
contour(xg,yg,zg,-1:1/18:2,'k') %corbes de nivell amb equidistancia 1/18
plot(0,0,'r*') %marquem els punts equilibri
plot(1,0,'r*')
plot([0,1],[0,0],'g') %segment que els uneix
axis equal
Totes les solucions amb condició inicial
essent 
són òrbites periòdiques i es recorren en sentit horari. Els seus períodes 
depenen de la condició inicial 
i varien entre 
(quan 
) i ∞ (quan 
).
són òrbites periòdiques i es recorren en sentit horari. Els seus períodes depenen de la condició inicial i varien entre (quan ) i ∞ (quan ).Càlcul d'un període amb l'opció Events
Ens proposem calcular, per a una condició inicial concreta, per exemple 
el període 
La idea és resoldre el PVI usant ode45 en un interval prou llarg 
de manera que 
i, per a la solució 
obtinguda, plantejar l'equació 
i buscar-ne la primera solució 
en què 
passi de positiu a negatiu.
el període La idea és resoldre el PVI usant ode45 en un interval prou llarg de manera que i, per a la solució obtinguda, plantejar l'equació i buscar-ne la primera solució en què passi de positiu a negatiu.Podem trobar 
directament usant ode45 amb l'opció Events, en la qual especificarem quina magnitud ha de canviar de signe (en aquest cas la corrdenada y) i en quin sentit (de positiu a negatiu o al revés). Ho farem amb una funció que anomenem, per exemple, condicio(t,X) i que podem definir al final del codi o bé en un fitxer apart (que s'ha d'anomenar condicio.m). D'aquesta manera obtenim, a més del vector t i la matriu X habituals, un vector te i una matriu Xe que ens donen, respectivament, els successius valors de la variable t per a les quals la condició es compleix, i els punts 
que corresponen als valors de t trobats. D'aquests, el període vindrà donat pel primer valor 
directament usant ode45 amb l'opció Events, en la qual especificarem quina magnitud ha de canviar de signe (en aquest cas la corrdenada y) i en quin sentit (de positiu a negatiu o al revés). Ho farem amb una funció que anomenem, per exemple, condicio(t,X) i que podem definir al final del codi o bé en un fitxer apart (que s'ha d'anomenar condicio.m). D'aquesta manera obtenim, a més del vector t i la matriu X habituals, un vector te i una matriu Xe que ens donen, respectivament, els successius valors de la variable t per a les quals la condició es compleix, i els punts que corresponen als valors de t trobats. D'aquests, el període vindrà donat pel primer valor opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicio);
%afegim, a les opcions habituals, l'opcio "Events" definida
%per una funcio que hem anomenat "condicio" (vegeu-la al final),
%en la qual demanem que y(t) canvii de signe, de positiu a negatiu
t0=0; X0=[0.5,0];
tf=25; %temps final (cal augmentar-lo si no es prou gran)
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
ne=length(te);
for i=1:ne
fprintf('te=%.6f\n',te(i)) %valors de t que compleixen la condicio
end
fprintf('\nperiode=%.6f\n',te(2)) %primera t>0 que compleix la condicio
Funció auxiliar
Podem definir la funció condicio al final del fitxer que conté el nostre codi, o bé en un fitxer apart (que s'haurà d'anomenar condicio.m).
function [expr,aturar,signe]=condicio(t,X)
expr=X(2);
%expressio de la qual volem detectar un canvi de signe,
%en aquest cas la y (segona component de X)
aturar=0;
%valor 1 si volem nomes el primer canvi de signe a partir de t0
%valor 0 si volem tots els canvis de signe entre t0 i tf
signe=-1;
%valor 1 si volem un canvi de signe de negatiu a positiu
%valor -1 si volem un canvi de signe de positiu a negatiu
%valor 0 si volem un canvi de signe qualsevol
end
________
(c) Numerical Factory (by Pere Gutiérrez)