Exercicis pràctica: EVOLUCIÓ D'UNA ÀREA EN UN SISTEMA 2D

Exercicis resolts

(A1) Considereu el sistema següent, que correspon a un pèndol amb fricció,
Prenem com a condicions inicials tots els punts de l'el·lipse de centre i semieixos parametritzada per
Integrant amb ode45 des de t0=0 fins a tf=5, calculeu per quin factor s'ha multiplicat la longitud de la corba de punts finals respecte la longitud de l'el·lipse inicial, usant la regla dels trapezis amb N=100 subintervals per a calcular les longituds. [ Nota: Com que el camp vectorial té divergència negativa, l'àrea encerclada per una corba tancada disminueix, però això pot no ser cert per a la longitud. ]
clear all
P=@(x,y) y;
Q=@(x,y) -sin(x)-0.1*y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
c=[0,2]; a=0.1; b=0.2;
N=100;
h=2*pi/N;
s=(0:h:2*pi)';
X0=[c(1)+a*cos(s),c(2)+b*sin(s)];
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
tf=5;
Xf=zeros(N+1,2);
for i=1:N+1
[t,X]=ode45(F,[0,tf],X0(i,:),opcions);
Xf(i,:)=X(end,:);
end
DX0=[-a*sin(s),b*cos(s)];
normaDX0=sqrt(DX0(:,1).^2+DX0(:,2).^2);
L0=trapz(s,normaDX0);
DXf=zeros(N+1,2);
for i=2:N
DXf(i,:)=(Xf(i+1,:)-Xf(i-1,:))/(2*h);
end
DXf(1,:)=(Xf(2,:)-Xf(N,:))/(2*h);
DXf(N+1,:)=DXf(1,:);
normaDXf=sqrt(DXf(:,1).^2+DXf(:,2).^2);
Lf=trapz(s,normaDXf);
factor=Lf/L0;
fprintf('factor = %.6f\n',factor)
factor = 10.261343
(A2) Considereu el camp vectorial i el seu sistema d'EDOs associat. Calculeu la circulació de al llarg de la trajectòria del sistema que compleix la condició inicial aplicant el mètode d'Euler i després la regla dels trapezis, les dues coses amb N=100 subintervals. [ Ind.: En aquest cas no cal usar derivació numèrica per a calcular ja que aquestes derivades vénen donades pel propi sistema d'EDOs. ]
clear all
P=@(x,y) -3*x+exp(x-y.^2);
Q=@(x,y) y+cos(x);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
t0=0; tf=4;
X0=[1,0];
N=100;
h=(tf-t0)/N;
[t,X]=odeEuler(F,[t0,tf],X0,h);
for i=1:N+1
DX(i,:)=F(t(i),X(i,:));
g(i)=dot(F(t(i),X(i,:)),DX(i,:));
end
circulacio=trapz(t,g);
fprintf('circulacio = %.6f\n',circulacio)
circulacio = 787.196264
(A3) Considerem el sistema d'EDOs
i prenem com a condicions inicials qualsevol punt de la circumferència Usant ode45, integrem el sistema d'EDOs des de t0=0 fins a tf=1.5. Calculeu la mitjana de la coordenada y sobre la corba dels punts finals, si apliquem la regla dels trapezis amb N=400 subintervals.
clear all
P=@(x,y) 1-y^2;
Q=@(x,y) x.^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
N=400;
h=2*pi/N;
s=(0:h:2*pi)';
X0=[cos(s),sin(s)];
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
tf=1.5;
Xf=zeros(N+1,2);
for i=1:N+1
[t,X]=ode45(F,[0,tf],X0(i,:),opcions);
Xf(i,:)=X(end,:);
end
DXf=zeros(N+1,2);
for i=2:N
DXf(i,:)=(Xf(i+1,:)-Xf(i-1,:))/(2*h);
end
DXf(1,:)=(Xf(2,:)-Xf(N,:))/(2*h);
DXf(N+1,:)=DXf(1,:);
normaDXf=sqrt(DXf(:,1).^2+DXf(:,2).^2);
Lf=trapz(s,normaDXf);
yf=Xf(:,2);
g=yf.*normaDXf;
intyf=trapz(s,g);
mitjana=intyf/Lf;
fprintf('mitjana = %.6f\n',mitjana)
mitjana = 0.633803
(A4) Donat el sistema d'EDOs
imposem condicions inicials prenent valors amb un increment h=0.005, i integrem fins a tf=4. Trobeu s per tal que el valor de sigui mínim.
clear all
P=@(x,y) -x+sin(y);
Q=@(x,y) y-cos(x.^2);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
h=0.005;
N=2/h;
s=(-1:h:1)';
X0=[s,0*s];
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
tf=4;
Xf=zeros(N+1,2);
for i=1:N+1
[t,X]=ode45(F,[0,tf],X0(i,:),opcions);
Xf(i,:)=X(end,:);
end
[ymin,i]=min(Xf(:,2));
smin=s(i);
fprintf('smin = %.2f\n',smin)
smin = 0.29
(A5) Considereu el sistema (lineal no homogeni a coeficients variables), donat per
Prenem condicions inicials al triangle de vèrtexs Integrant des de t0=0 fins a tf=2, calculeu per quin factor s'ha multiplicat l'àrea del triangle final respecte l'àrea del triangle inicial. [ Ind.: Com que en un sistema lineal les rectes es transformen en rectes, i per tant els triangles en triangles, només cal trobar els punts finals per als 3 vèrtexs del triangle, i aplicar la fórmula de l'àrea del triangle. ]
[ Comprovació: Per la fórmula de Liouville, el quocient entre les dues àrees és ]
A=@(t) [sin(t),t/2; 3*exp(-t),0];
b=@(t) [-t.*cos(t); 1];
F=@(t,X) A(t)*X+b(t);
%definicio alternativa:
% P=@(t,x,y) sin(t).*x+t/2.*y-t.*cos(t);
% Q=@(t,x,y) 3*exp(-t).*x+1;
% F=@(t,X) [P(t,X(1),X(2));Q(t,X(1),X(2))];
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
tf=2;
X0=[1,0; 0,1; -0.5,0];
Xf=zeros(3,2);
for i=1:3
[t,X]=ode45(F,[0,tf],X0(i,:),opcions);
Xf(i,:)=X(end,:);
end
area0=1/2*det([X0(2,:)-X0(1,:);X0(3,:)-X0(1,:)]);
areaf=1/2*det([Xf(2,:)-Xf(1,:);Xf(3,:)-Xf(1,:)]);
factor=areaf/area0;
fprintf('factor = %.6f\n',factor)
factor = 4.121210
% Comprovacio:
factorexacte=exp(-cos(tf)+1); %formula de Liouville
dif=factor-factorexacte;
fprintf('\ncomprovacio: dif = %.4e\n',dif)
comprovacio: dif = 7.5290e-09

Altres exercicis

(B1) Sigui Cf(s)=(thetaf(s),omegaf(s)), 1<=s<=4, l'evolució del segment C0(s)=(theta0(s),omega0(s))=(0,s) pel sistema del pèndol sense fricció theta' = omega, omega' = -alpha sin(theta), durant l'interval de temps [t0,tf]=[0,5] amb alpha=2.00. Calculeu Cf(sv), on sv=1:0.01:4, usant la instrucció ode45. Després, trobeu el paràmetre s que correspon al punt de Cf(sv) més proper al punt (pi,0) (és a dir, a la posició d'equilibri superior). [ Comprovació: Avaluant l'energia E = omega^2/2 - alpha cos(theta) sobre les components del punt obtingut Cf(s)=(theta,omega), tindrem E ~ alpha, que és l'energia del punt (pi,0). ]
[ Solució: 2.83 ]
(B2) Considereu el camp vectorial F(x,y,z) = (a x + y^2, x-y+z, y+exp(-z)) i el seu sistema d'EDOs associat. Calculeu la circulació de F al llarg de la trajectòria del sistema (x(t),y(t),z(t)), amb t entre 0 i b, que compleix la condició inicial (x(0),y(0),z(0))=(1,0,-1), aplicant el mètode d'Euler i després la regla dels trapezis, ambdós amb N=50 intervals. Preneu a=1.50, b=1.00.
[ Solució: 30.435188 ]
(B3) Sigui (x(t),y(t)) la solució periòdica del sistema depredador-presa x'=x(1-y), y'=y(x-1), tal que x(0)=x0 i y(0)=1, essent x0=3.00. Volem calcular una aproximació del període T, del qual sabem que 2pi<T<4pi. Calculeu la solució amb la instrucció ode45 (amb tolerància 1e-8 per a l'error absolut i relatiu) durant un interval de temps [t0,tf] = [0,4pi]. La instrucció ode45 ens dóna un vector columna t amb tots els instants on tenim la solució aproximada. Busqueu l'instant tm, amb 2 pi < tm < 4 pi, per al qual la solució és més a prop de la posició inicial (x0,1). Llavors tenim T ~ tm.
[ Solució: 7.261048 ]
(B4) Considereu el sistema x' = y^2 - x^2, y' = d x + y - 2. Considereu l'el·lipse de centre (-2,1.5) i semieixos a i b, parametritzada per (-2 + a cos(s), 1.5 + b sin(s)), amb s entre 0 i 2pi. Prenem N=100 punts equiespaiats en s sobre l'el·lipse com a condicions inicials i integrem el sistema per a temps entre t0=0 i tf=0.25. Tenim l'àrea inicial A0 limitada per l'el·lipse, i la final Af limitada pels punts a temps tf. Calculeu el quocient Af/A0. Preneu d=1.30, a=0.20, b=0.40.
[ Solució: 4.396829 ]
(B5) Considereu el sistema d'EDOs x' = x+y, y' = 2x + c sin(y) i preneu com a condicions inicials (x(0),y(0)) qualsevol punt de l'el·lipse de centre (1,2) i semieixos a, b, amb parametrització (1 + a cos(s), 2 + b sin(s)), amb 0<=s<=2pi. Usant ode45 integreu el sistema d'EDOs des de t0=0 fins a tf. Calculeu la suma de la mitjana de les coordenades x i y sobre la corba de punts finals, si apliquem la regla dels trapezis amb N=300 intervals. Preneu c=0.50, a=1.00, b=2.00, tf=0.50.
[ Solució: 7.443981e+00 ]
(B6) Considereu el sistema d'EDOs x' = a x + y, y' = x-y^2 i preneu com a condicions inicials (x(0),y(0)) qualsevol punt del tros de recta que uneix els punts (1,0) i (2,2), parametritzada per (s,2s-2), amb s entre 1 i 2. Usant ode45 integreu el sistema d'EDOs des de t0=0 fins a tf. Calculeu la mitjana de la coordenada x sobre la corba de punts finals, si apliquem la regla dels trapezis amb N=200 intervals. Preneu a=0.80, tf=1.10. [ Ind.: Com que no tenim una corba tancada, per poder calcular numèricament la derivada als extrems haurem de considerar valors de s entre 1-h i 2+h, essent h=1/N. ]
[ Solució: 5.828640 ]
(B7) Considereu el sistema d'EDOs x'=x+y, y'=x+ky+5 i preneu com a condicions inicials qualsevol punt de la circumferència de centre (a,b) i radi r, parametritzada per (a + r cos(s), b +r sin(s)), amb s entre 0 i 2pi. Usant ode45, integreu el sistema des de t0=0 fins a tf=0.5, obtenint així una nova corba. Aixequem, sobre la nova corba, una tanca d'altura definida per la funció f(x,y)=x+y. Aplicant trapezis amb N=100 intervals, calculeu-ne l'àrea. Preneu els valors k=1.30, a=3, b=0, r=1.
[ Solució: 173.394584 ]
(B8) Sigui (theta(t),omega(t)) la solució periòdica del pèndol sense fricció theta' = omega, omega' = -sin(theta), tal que theta(0) = a pi i omega(0) = 0, essent a=0.50. Volem calcular una aproximació del període T, del qual sabem que 2pi<T<4pi. Per la simetria de les solucions, és suficient buscar el semiperíode T/2. Calculeu la solució amb la instrucció ode45 (amb tolerància 1e-8 per a l'error absolut i relatiu) durant un interval de temps [t0,tf] = [0,2pi]. La instrucció ode45 ens dóna un vector columna t amb tots els instants on tenim la solució aproximada. Busqueu l'instant tm, amb t0<tm<tf, en què el pèndol és més a prop de l'angle oposat a l'inicial, és a dir, busqueu tm tal que theta(tm) ~ -theta(0). Llavors el període és T ~ 2 tm. [ Comprovació: (theta(T),omega(T)) ~ (theta(0),0). ]
[ Solució: 7.415136 ]
(B9) Donat el sistema d'EDOs x'=-y+sin(x), y'=x+cos(y), preneu com a condicions inicials punts sobre el segment parametritzat per (s,0) amb s pertanyent a l'interval [1,2] amb increment h=0.01. Utilitzeu ode45 per obtenir numèricament la solució per a cada s, per a temps entre t0=0 i tf=3. Trobeu el màxim dels valors de x(3) obtinguts d'aquesta manera.
[ Solució: -2.483834 (obtingut per a s=1.51) ]