Pràctica: DETERMINACIÓ D'ESDEVENIMENTS EN EDOs (versió ampliada)
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.
En primer lloc, veurem un mètode genèric de càlcul de zeros, com és el mètode de Newton, i com aplicar-lo a la solució d'una EDO obtinguda amb ode45. Després veurem com obtenir-ho directament utilitzant ode45 amb l'opció Events.
Càlcul de zeros: la funció fzero i el mètode de Newton
Donada una funció d'una variable 
volem trobar el seus zeros, és a dir resoldre l'equació 
Abans de tot, esmentem la funció de Matlab fzero, que incorpora un mètode numèric (mètode de Dekker-Brent, una combinació de bisecció i secant), i que podem usar de dues maneres diferents:
volem trobar el seus zeros, és a dir resoldre l'equació Abans de tot, esmentem la funció de Matlab fzero, que incorpora un mètode numèric (mètode de Dekker-Brent, una combinació de bisecció i secant), i que podem usar de dues maneres diferents:- fzero(g,tini), si coneixem una aproximació inicial
del zero; - fzero(g,[a,b]), si coneixem un interval
en què es produeix un canvi de signe.
Estem suposant que la funció ve definida per una ordre del tipus g=@(t).... En el cas d'una funció definida a través d'un fitxer gfun.m, caldria canviar g per @gfun.
D'aquesta manera només podem determinar un zero. Per trobar-los tots, cal fer un escombrat previ de valors de la funció per tal de detectar canvis de signe (amb un increment suficientment petit de la variable 
per no saltar-nos cap zero).
per no saltar-nos cap zero).Per exemple, la funció 
té 3 zeros, un dels quals és proper a 
i pertany a l'interval 
en què es produeix un canvi de signe. Representem la funció i vegem com trobar aquest zero amb fzero, de les dues maneres.
té 3 zeros, un dels quals és proper a i pertany a l'interval en què es produeix un canvi de signe. Representem la funció i vegem com trobar aquest zero amb fzero, de les dues maneres.g=@(t) t+4*cos(t);
t=-4:0.01:4;
plot(t,g(t))
grid
z1=fzero(g,2.1)
z1 = 2.1333
z2=fzero(g,[0,2.5])
z2 = 2.1333
Vegem ara el mètode de Newton. Per trobar un zero α de la funció 
partint d'una aproximació inicial 
la idea del mètode és construir aproximacions successives 
i cada nova aproximació 
s'obté de l'anterior 
com a intersecció de la recta tangent a la gràfica en el punt 
amb l'eix d'abscisses. Això ens dóna la fórmula iterativa següent,
partint d'una aproximació inicial la idea del mètode és construir aproximacions successives i cada nova aproximació s'obté de l'anterior com a intersecció de la recta tangent a la gràfica en el punt amb l'eix d'abscisses. Això ens dóna la fórmula iterativa següent,
En general, les aproximacions obtingudes convergiran ràpidament al zero α si escollim l'aproximació inicial 
prou propera al zero. De fet, si 
(zero simple) l'ordre d'aproximació és quadràtic: escrivint 
es compleix que 
Això ens ve a dir que a cada iteració doblem el nombre de xifres decimals correctes. En canvi, si 
(zero múltiple) l'aproximació no és quadràtica sinó lineal.
prou propera al zero. De fet, si (zero simple) l'ordre d'aproximació és quadràtic: escrivint es compleix que Això ens ve a dir que a cada iteració doblem el nombre de xifres decimals correctes. En canvi, si (zero múltiple) l'aproximació no és quadràtica sinó lineal.Per aturar les iteracions, podem escollir com a criteri que la diferència 
entre les dues darreres aproximacions calculades sigui més petita que una tolerància d'error tol prefixada.
entre les dues darreres aproximacions calculades sigui més petita que una tolerància d'error tol prefixada.Com a exemple, vegem com calcular un dels zeros de la funció 
a partir de l'aproximació inicial 
amb una tolerància d'error tol=1e-12.
a partir de l'aproximació inicial amb una tolerància d'error tol=1e-12.g=@(t) t+4*cos(t);
Dg=@(t) 1-4*sin(t); %funcio i derivada
tini=2.1; %aproximacio inicial
fprintf('n=0, t=%.15f\n',tini)
tol=1e-12;
nmax=20; %nombre maxim d'iteracions permeses
n=0; %comptador d'iteracions
incr=1; %valor fals per obligar a fer la primera iteracio
tprev=tini;
while( abs(incr)>=tol && n<=nmax )
n=n+1;
incr=-g(tprev)/Dg(tprev);
tnou=tprev+incr; %nova aproximacio
fprintf('n=%d, t=%.15f, incr=%.6e\n',n,tnou,incr)
tprev=tnou; %redefinim per a la iteracio seguent
end
Aplicació: període d'una òrbita periòdica d'un sistema d'EDOs
Considerem el sistema d'EDOs
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).
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).close all %tanquem figures anteriors (optatiu)
figure %nova figura
hold on
axis equal
V=@(x,y) y.^2/2+x.^2/2-x.^3/3;
x=-0.7:0.01:1.5; %finestra a visualitzar
y=-0.8:0.01:0.8;
[X,Y]=meshgrid(x,y);
Z=V(X,Y);
contour(X,Y,Z,-1:0.05:2,'k') %corbes de nivell amb equidistancia 0.05
plot(0,0,'r*') %marquem els punts equilibri
plot(1,0,'r*')
plot([0,1],[0,0],'g') %segment que els uneix
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 ).Ens proposem calcular, per a una condició inicial concreta, per exemple 
el període 
Primer resoldrem el PVI usant ode45 en un interval prou llarg 
de manera que 
i, per a la solució 
obtinguda, plantejarem l'equació
el període Primer resoldrem el PVI usant ode45 en un interval prou llarg de manera que i, per a la solució obtinguda, plantejarem l'equació
i buscarem la primera solució 
en què 
passi de positiu a negatiu.
en què passi de positiu a negatiu.En aquest cas, no podem usar fzero ja que no coneixem la fórmula explícita de 
però sí podem aplicar el mètode de Newton, calculant una aproximació de amb ode45 i usant que 
Trobarem l'aproximació inicial per iniciar les iteracions del mètode de Newton a partir de la d'una primera taula de valors obtinguda amb ode45.
però sí podem aplicar el mètode de Newton, calculant una aproximació de amb ode45 i usant que Trobarem l'aproximació inicial per iniciar les iteracions del mètode de Newton a partir de la d'una primera taula de valors obtinguda amb ode45.opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
P=@(x,y) y;
Q=@(x,y) -x+x.^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
t0=0; X0=[0.5,0];
tf=25; %temps final (cal augmentar-lo si no es prou gran)
[t,X]=ode45(F,[t0,tf],X0,opcions);
nt=length(t); %nombre de punts obtinguts
% Dins de la segona columna de la matriu X obtinguda (que ens dona y(t)),
% busquem el primer canvi de signe de positiu a negatiu
i=1;
while( ~(X(i,2)>0 && X(i+1,2)<0) ) %iterem mentre no es compleixi
i=i+1;
end
[yini,j]=min(X(i:i+1,2));
tini=t(i+j-1);
Xini=X(i+j-1,:);
%entre els valors t(i),t(i+1) hem escollit
%el que correspon a abs(y(t)) minima
fprintf('n=0, t=%.8f\n',tini)
% Apliquem el metode de Newton
tol=1e-6;
%no cal una tolerancia molt petita, ja que l'EDO no es resol exactament
nmax=20;
n=0;
incr=1; %valor fals
tprev=tini;
Xprev=Xini;
DXprev=F(tprev,Xprev);
while( abs(incr)>=tol && n<=nmax )
n=n+1;
incr=-Xprev(2)/DXprev(2); %segones components: y(t), y'(t)
tnou=tprev+incr; %nova aproximacio
fprintf('n=%d, t=%.8f, incr=%.6e\n',n,tnou,incr)
[t,X]=ode45(F,[tprev,tnou],Xprev,opcions);
tprev=tnou; %redefinim per a la iteracio seguent
Xprev=X(end,:);
DXprev=F(tprev,Xprev);
end
fprintf('\nperiode=%.6f\n',tnou)
Obtenció del període amb l'opció Events
Una manera més de trobar el període a l'exemple anterior, és usar ode45 amb l'opció Events, en la qual especificarem quina magnitud ha de canviar de signe i en quin sentit (negatiu a positiu 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 
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
P=@(x,y) y;
Q=@(x,y) -x+x.^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
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 definirla 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)