ALGUNES QÜESTIONS D'EXÀMENS ANTICS DE MATLAB - Equacions Diferencials

(Maig 2024)

P1: Integració sobre corbes

[C11] Considerem la cardioide definida en coordenades polars per

amb densitat donada per

Trobeu la primera coordenada

del seu centre de masses

clear
 
r=@(t) 2.50*(1+cos(t));
Dr=@(t) -2.50*sin(t);
nt=@(t) sqrt(r(t).^2+Dr(t).^2);
 
x=@(t) r(t).*cos(t);
y=@(t) r(t).*sin(t);
rho=@(y) 0.60+y.^2;
 
g=@(t) rho(y(t)).*nt(t);
massa=integral(g,0,2*pi);
h=@(t) x(t).*rho(y(t)).*nt(t);
intx=integral(h,0,2*pi);
xc=intx/massa;
 
fprintf('xCM = %.6f\n',xc)
xCM = 1.962475

[C12] Calculeu la circulació del camp

essent

al llarg de la corba C (la vora d'un pètal de rosa), definida en coordenades polars

per l'equació

i recorreguda en sentit antihorari, usant quadratura adaptativa.

clear
 
beta=2.60;
P=@(x,y) beta*cos(x.*y);
Q=@(x,y) log(beta+x+y);
 
r=@(t) 1.50*sin(3*t);
x=@(t) r(t).*cos(t);
y=@(t) r(t).*sin(t);
 
Dr=@(t) 4.50*cos(3*t);
Dx=@(t) Dr(t).*cos(t)-r(t).*sin(t);
Dy=@(t) Dr(t).*sin(t)+r(t).*cos(t);
 
g=@(t) P(x(t),y(t)).*Dx(t)+Q(x(t),y(t)).*Dy(t);
circulacio=integral(g,0,pi/3);
 
fprintf('circulacio = %.6f\n',circulacio)
circulacio = 0.624416

[C13] Calculeu la massa del filferro

si la densitat a cada punt ve donada pel cub de la distància a l'eix x (usant quadratura adaptativa).

clear
 
y=@(x) x.^2-3.50;
Dy=@(x) 2*x;
nn=@(x) sqrt(1+Dy(x).^2);
rho=@(x,y) abs(y).^3;
 
g=@(x) rho(x,y(x)).*nn(x);
massa=integral(g,0,3.90);
 
fprintf('massa = %.6e\n',massa)
massa = 4.803010e+03

[C14] Calculeu la circulació del camp

al llarg de la corba orientada parametritzada per

usant quadratura adaptativa.

clear
 
P=@(x,y) 2.50*x.^3-2.97*y;
Q=@(x,y) 1.03*y.*x.^2;
 
x=@(t) 1+cos(exp(t.^2));
y=@(t) t.^3.*sin(3.80*t);
Dx=@(t) -2*sin(exp(t.^2)).*exp(t.^2).*t;
Dy=@(t) 3*t.^2.*sin(3.80*t)+3.80*t.^3.*cos(3.80*t);
 
g=@(t) P(x(t),y(t)).*Dx(t)+Q(x(t),y(t)).*Dy(t);
circulacio=integral(g,0,pi);
 
fprintf('circulacio = %.6f\n',circulacio)
circulacio = 315.826406

[C15] Considerem la corba parametritzada

Calculeu la circulació del camp

usant quadratura adaptativa.

clear
 
P=@(x,y,z) x-0.90*z.^2;
Q=@(x,y,z) 0.70*y;
R=@(x,y,z) 0.90*z;
 
x=@(t) exp(1.20*t).*cos(t);
y=@(t) exp(-1.20*t).*sin(t);
z=@(t) t.^2;
 
Dx=@(t) exp(1.20*t).*(1.20*cos(t)-sin(t));
Dy=@(t) exp(-1.20*t).*(cos(t)-1.20*sin(t));
Dz=@(t) 2*t;
 
g=@(t) P(x(t),y(t),z(t)).*Dx(t)+Q(x(t),y(t),z(t)).*Dy(t)...
    +R(x(t),y(t),z(t)).*Dz(t);
circul=integral(g,0,2*pi);
 
fprintf('circulacio = %.6e\n',circul)
circulacio = -1.279589e+05

[C16] Considerem la corba C del pla, definida en coordenades polars per

és el centre de masses de

calculeu el moment d'inèrcia de C respecte la recta

suposant densitat constant

Useu quadratura adaptativa per fer els càlculs.

clear
 
r=@(t) 2.83*sin(3*t)+0.05*sin(30*t).^2;
Dr=@(t) 2.83*3*cos(3*t)+0.05*60*sin(30*t).*cos(30*t);
x=@(t) r(t).*cos(t);
%y=@(t) r(t).*sin(t);
Dx=@(t) Dr(t).*cos(t)-r(t).*sin(t);
Dy=@(t) Dr(t).*sin(t)+r(t).*cos(t);
 
t0=0; tf=pi/3;
 
nt=@(t) sqrt(Dx(t).^2+Dy(t).^2);
long=integral(nt,t0,tf);
 
h1=@(t) x(t).*nt(t);
intx=integral(h1,t0,tf);
x0=intx/long; %coordenada x del centre
 
rho=1.22;
h2=@(t) (x(t)-x0).^2.*nt(t); % |x(t)-x0| es la distancia a la recta x=x0
I=rho*integral(h2,t0,tf);
 
fprintf('I = %.6f\n',I)
I = 4.814698

[C17] Calculeu la circulació del camp

al llarg de l'el·lipse

orientada en sentit antihorari, essent

usant quadratura adaptativa.

clear
 
a=0.60; b=2.20;
 
P=@(x,y) exp(x).*y.^2;
Q=@(x,y) cos(x.*y+2*y);
 
x=@(t) a*cos(t); y=@(t) b*sin(t);
Dx=@(t) -a*sin(t); Dy=@(t) b*cos(t);
 
g=@(t) P(x(t),y(t)).*Dx(t)+Q(x(t),y(t)).*Dy(t);
circulacio=integral(g,0,2*pi);
 
fprintf('circulacio = %.6f\n',circulacio)
circulacio = -0.998390

[C18] Donada la corba del pla

amb densitat de massa

calculeu el seu moment d'inèrcia al voltant de l'eix x usant quadratura adaptativa. Preneu els valors

clear
 
a=1.30; b=6.00;
y=@(x) a+x.^2.*sin(x);
Dy=@(x) 2*x.*sin(x)+x.^2.*cos(x);
rho=@(x,y) x.^2+3*y+b;
 
nt=@(x) sqrt(1+Dy(x).^2);
g=@(x) y(x).^2.*rho(x,y(x)).*nt(x);
I=integral(g,0,pi);
 
fprintf('I = %.6e\n',I)
I = 2.454626e+03

[C19] Calculeu la mitjana de la funció

essent

al llarg de la semiel·lipse parametritzada per

amb

aplicant la regla dels trapezis amb

subintervals.

clear
 
a=2.50; b=1.70;
f=@(x,y) 1./(a*x.^2+y.^4);
x=@(t) b*cos(t); y=@(t) sin(t);
Dx=@(t) -b*sin(t); Dy=@(t) cos(t);
 
N=1000; h=pi/N; tv=0:h:pi;
xv=x(tv); yv=y(tv);
Dxv=Dx(tv); Dyv=Dy(tv);
 
normav=sqrt(Dxv.^2+Dyv.^2);
longaprox=trapz(tv,normav);
 
fv=f(xv,yv);
gv=fv.*normav;
intaprox=trapz(tv,gv);
 
mitjanaaprox=intaprox/longaprox;
fprintf('mitjanaaprox = %.6f\n',mitjanaaprox)
mitjanaaprox = 0.430075

P2: Integració sobre superfícies

[S06] Considerem la superfície d'equació

i tallada pel cilindre

sobre la qual la temperatura ve donada per

Calculeu-ne la temperatura mitjana.

clear
 
z=@(x,y) x.^2-y.^2;
Dzx=@(x,y) 2*x;
Dzy=@(x,y) -2*y;
 
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
nnpolars=@(r,th) nn(r.*cos(th),r.*sin(th)).*r;
massa=integral2(nnpolars,0,1,0,2*pi);
 
T=@(x,y,z) 2.10*x.^2+1.70*y.^2+3.10*z;
g=@(x,y) T(x,y,z(x,y)).*nn(x,y);
gpolars=@(r,th) g(r.*cos(th),r.*sin(th)).*r;
intg=integral2(gpolars,0,1,0,2*pi);
 
Tm=intg/massa;
fprintf('Tm = %.6f\n',Tm)
Tm = 1.061981

[S07] Calculeu la massa de la superfície definida per

si la densitat ve donada per

clear
 
z=@(x,y) 3.40-x.^2-y.^2;
Dzx=@(x,y) -2*x;
Dzy=@(x,y) -2*y;
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
rho=@(x,y,z) (x-1).^2 + y.*z;
 
g=@(x,y) rho(x,y,z(x,y)).*nn(x,y);
gpolars=@(r,th) g(r.*cos(th),r.*sin(th)).*r;
 
rmin=sqrt(3.40-0.60); rmax=sqrt(3.40);
  %corona circular de radis rmin i rmax
massa=integral2(gpolars,rmin,rmax,0,2*pi);
fprintf('massa = %.6f\n',massa)
massa = 17.604748

[S08] Calculeu la massa de la superfície

tallada pel cilindre

si la densitat superficial és proporcional al quadrat de la distància al pla

essent

la constant de proporcionalitat.

clear
 
h=@(x,y) x.^3-y.^3;
Dhx=@(x,y) 3.*x.^2; Dhy=@(x,y) -3.*y.^2;
nn=@(x,y) sqrt(1+Dhx(x,y).^2+Dhy(x,y).^2);
rho=@(x,y,z) 4.70*z.^2;
 
g=@(x,y) rho(x,y,h(x,y)).*nn(x,y); %funcio a integrar
gpolars=@(r,th) g(r.*cos(th),r.*sin(th)).*r;
massa=integral2(gpolars,0,sqrt(1.44),0,2*pi);
 
fprintf('massa = %.6f\n',massa)
massa = 32.626143

[S09] Trobeu la tercera coordenada

del centre geomètric de la superfície

clear
 
z=@(x,y) 1+1.90*x.^4+2.80*y.^2;
Dzx=@(x,y) 4*1.90*x.^3;
Dzy=@(x,y) 2*2.80*y;
 
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
area=integral2(nn,-1,1,-1,1);
 
g=@(x,y) z(x,y).*nn(x,y);
I=integral2(g,-1,1,-1,1);
 
zc=I/area;
fprintf('zCG = %.6f\n',zc)
zCG = 2.735052

[S10] Trobeu la massa de la porció de paraboloide

amb densitat

clear
 
a2=2.89; b2=6.25;
z=@(x,y) x.^2/a2+y.^2/b2;
Dzx=@(x,y) 2*x/a2; Dzy=@(x,y) 2*y/b2;
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
rho=@(x,y,z) 1+x.^2.*z;
 
g=@(x,y) rho(x,y,z(x,y)).*nn(x,y);
a=sqrt(a2); b=sqrt(b2); %semieixos ellipse
gpolars=@(r,th) g(a*r.*cos(th),b*r.*sin(th)).*(a*b*r);
  %hem passat a polars adaptades i multiplicat pel jacobia
massa=integral2(gpolars,0,1,0,2*pi);
 
fprintf('massa = %.6f\n',massa)
massa = 25.126766

[S11] Calculeu el flux del camp

a través de la superfície

essent

orientada per la normal exterior.

clear
 
P=@(x,y,z) y.*z;
Q=@(x,y,z) 1.60*x.^4.*exp(x.^2);
R=@(x,y,z) y.^2;
 
b=1.10;
 
% aplicant el teorema de Gauss:
%   tanquem S per un tros de pla (ellipse inclinada)
%   que es projecta sobre el cercle de centre (b/2,0) i radi b/2;
%   com que div(F)=0 nomes hem de calcular el flux sobre el tros de pla
%   (orientat cap amunt)
 
h=@(x,y) b-x; %el pla vist com a grafica
g1=@(x,y) -P(x,y,h(x,y))+R(x,y,h(x,y));
%producte escalar de F amb el vector normal del pla (-1,0,1)
 
g1polars=@(r,th) g1(b/2+r.*cos(th),r.*sin(th)).*r;
flux1=integral2(g1polars,0,b/2,0,2*pi);
fprintf('flux1 = %.6f\n',flux1)
flux1 = 0.071869
 
% sense aplicar Gauss:
%   parametritzant el tros d'ellipsoide com a grafica,
%   amb domini el mateix cercle de centre (b/2,0) i radi b/2
 
z=@(x,y) sqrt(b^2-x.^2+2*y.^2);
Dzx=@(x,y) -x./z(x,y);
Dzy=@(x,y) -2*y./z(x,y);
g2=@(x,y) -P(x,y,z(x,y)).*Dzx(x,y)-Q(x,y,z(x,y)).*Dzy(x,y)+R(x,y,z(x,y));
 
g2polars=@(r,th) g2(b/2+r.*cos(th),r.*sin(th)).*r;
flux2=integral2(g2polars,0,b/2,0,2*pi);
fprintf('flux2 = %.6f\n',flux2)
flux2 = 0.071869

[S12] Donat el camp vectorial

essent

i considerant l'esfera centrada a l'origen i de radi

trobeu el flux de

a través del casquet definit per la latitud

orientat pel vector normal exterior.

clear
 
beta=0.70;
P=@(x,y,z) x.*z.*cos(beta*y.^2);
Q=@(x,y,z) exp(beta*y.*z);
R=@(x,y,z) y.^2.*sin(beta*x.^2);
 
A=1.80; %radi
 
%parametritzant amb esferiques
x=@(th,phi) A*cos(phi).*cos(th);
y=@(th,phi) A*cos(phi).*sin(th);
z=@(th,phi) A*sin(phi);
 
nx=@(th,phi) A*cos(phi).*x(th,phi);
ny=@(th,phi) A*cos(phi).*y(th,phi);
nz=@(th,phi) A*cos(phi).*z(th,phi);
  %el vector normal es (nx,ny,xz)=A*cos(phi)*(x,y,z)
 
g1=@(th,phi)...
  P(x(th,phi),y(th,phi),z(th,phi)).*nx(th,phi)...
  +Q(x(th,phi),y(th,phi),z(th,phi)).*ny(th,phi)...
  +R(x(th,phi),y(th,phi),z(th,phi)).*nz(th,phi);
flux1=integral2(g1,0,2*pi,pi/4,pi/2);
fprintf('flux1 = %.6f\n',flux1)
flux1 = 4.000642
 
%parametritzant com a grafica
zgr=@(x,y) sqrt(A^2-x.^2-y.^2);
Dzx=@(x,y) -x./zgr(x,y);
Dzy=@(x,y) -y./zgr(x,y);
g2=@(x,y)...
  -P(x,y,zgr(x,y)).*Dzx(x,y)-Q(x,y,zgr(x,y)).*Dzy(x,y)+R(x,y,zgr(x,y));
g2polars=@(r,th) g2(r.*cos(th),r.*sin(th)).*r;
flux2=integral2(g2polars,0,A/sqrt(2),0,2*pi);
fprintf('flux2 = %.6f\n',flux2)
flux2 = 4.000642

[S13] Sigui S la superfície de revolució definida per l'equació

i limitada per les desigualtats

Trobeu la coordenada

del centre geomètric

clear
 
xc=@(z) sqrt(1+2.10*z.^2); %corba generatriu (tros d'hiperbola)
Dxc=@(z) 2.10*z./xc(z);
 
xs=@(th,z) xc(z).*cos(th); %sup. de revolucio (hiperboloide 1 fulla)
ys=@(th,z) xc(z).*sin(th);
 
nn=@(z) xc(z).*sqrt(Dxc(z).^2+1);
area=pi*integral(nn,-5.70,5.70); %meitat sup. de revolucio
 
ydS=@(th,z) ys(th,z).*nn(z);
inty=integral2(ydS,0,pi,-5.70,5.70);
 
yCG=inty/area;
fprintf('yCG = %.6f\n',yCG)
yCG = 3.540827

[S14] Calculeu la càrrega total de la superfície definida per

si la densitat de càrrega a cada punt és

clear
 
z=@(x,y) 4.20+x.^2+y.^2;
Dzx=@(x,y) 2*x; Dzy=@(x,y) 2*y;
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
rho=@(x,y,z) y.^2+z;
 
g=@(x,y) rho(x,y,z(x,y)).*nn(x,y);
gpolars=@(r,th) g(r.*cos(th),r.*sin(th)).*r;
carrega=integral2(gpolars,0,sqrt(7.70-4.20),0,2*pi);
 
fprintf('carrega = %.6f\n',carrega)
carrega = 216.891475

[S15] Calculeu la tercera coordenada del centre geomètric

de la superfície de revolució obtinguda fent girar la corba parametritzada per

al voltant de l'eix vertical

[ Ind.: recordeu que, en el cas d'una superfície de revolució, l'element de superfície ve donat per

]

clear
 
r=@(t) t;
z=@(t) sqrt(1+7.20*t);
Dz=@(t) 3.60./sqrt(1+7.20*t);
 
elem=@(t) sqrt(1+Dz(t).^2).*r(t);
Area=2*pi*integral(elem,0,2.60);
 
zelem=@(t) z(t).*elem(t);
zc=2*pi*integral(zelem,0,2.60)/Area;
 
fprintf('zCG = %.6f\n',zc)
zCG = 3.519907

[S16] Calculeu el flux del camp vectorial

a través de la superfície

pertany al triangle de vèrtexs

i la superfície ve orientada pel vector normal que té la tercera component positiva.

clear
 
P=@(x,y,z) sin(x.^2.*y);
Q=@(x,y,z) cos(x.*z);
R=@(x,y,z) x+y;
 
z=@(x,y) 4.10*x.^2-1.50*y.^2;
Dzx=@(x,y) 8.20*x;
Dzy=@(x,y) -3.00*y;
 
g=@(x,y) -P(x,y,z(x,y)).*Dzx(x,y)-Q(x,y,z(x,y)).*Dzy(x,y)+R(x,y,z(x,y));
ymin=@(x) -1.80*x;
ymax=@(x) 3.10*x;
flux=integral2(g,0,1,ymin,ymax);
 
fprintf('flux = %.6f\n',flux)
flux = 0.498837

[S17] Sigui S la porció del paraboloide

definida per

Calculeu la massa de

sabent que la seva densitat superficial és

clear
 
h=@(x,y) 1.20*x.^2+1.50*y.^2;
Dhx=@(x,y) 2.40*x;
Dhy=@(x,y) 3.00*y;
nn=@(x,y) sqrt(1+Dhx(x,y).^2+Dhy(x,y).^2);
 
rho=@(x,y,z) 2*z.^2;
 
g=@(x,y) rho(x,y,h(x,y)).*nn(x,y);
x2=1/0.90;
ymax=@(x) (1-0.90*x)/0.60; %funcio
massa=integral2(g,0,x2,0,ymax);
 
fprintf('massa = %.6f\n',massa)
massa = 9.400789

[S18] Sigui S la porció de l'el·lipsoide d'equació

delimitat per

Prenent els valors

calculeu el flux del camp vectorial

a través de

orientada pel vector normal

al punt

[ Ind.: hi ha diverses maneres de fer l'exercici; potser la més senzilla és escriure la superfície com una gràfica

]

clear
 
a=2.90; b=1.70;
h=@(x,y) sqrt(1-x.^2/a^2-y.^2/b^2);
Dhx=@(x,y) (-x/a^2)./h(x,y);
Dhy=@(x,y) (-y/b^2)./h(x,y);
P=@(x,y,z) sin(x.^2.*y);
Q=@(x,y,z) y;
R=@(x,y,z) z;
g=@(x,y) -P(x,y,h(x,y)).*Dhx(x,y)-Q(x,y,h(x,y)).*Dhy(x,y)+R(x,y,h(x,y));
 
% fent la integral doble iterada:
ymin=@(x) -b*sqrt(1-x.^2/a^2);
ymax=@(x) b*sqrt(1-x.^2/a^2);
flux1=integral2(g,-a,a,ymin,ymax);
fprintf('flux1 = %.6f\n',flux1)
flux1 = 20.650736
 
% fent la integral doble per polars adaptades:
gpolars=@(r,th) g(a*r.*cos(th),b*r.*sin(th)).*(a*b*r);
flux2=integral2(gpolars,0,1,0,2*pi);
fprintf('flux2 = %.6f\n',flux2)
flux2 = 20.650736
 
% usant el teorema de Gauss, i fent la integral triple
% per esferiques adaptades:
divF=@(x,y,z) 2*x.*cos(x.^2.*y)+2;
divFesf=@(r,th,phi)...
  divF(a*r.*cos(phi).*cos(th),b*r.*cos(phi).*sin(th),r.*sin(phi))...
  .*(a*b*r.^2.*cos(phi));
flux3=integral3(divFesf,0,1,0,2*pi,0,pi/2);
fprintf('flux3 = %.6f\n',flux3)
flux3 = 20.650736

[S19] Trobeu l'àrea de la porció del semiel·lipsoide

que es troba per sobre del triangle del pla

de vèrtexs

clear
 
arrel=@(x,y) sqrt(1-x.^2/4-y.^2/9);
Dzx=@(x,y) -1.10*x./(4*arrel(x,y));
Dzy=@(x,y) -1.10*y./(9*arrel(x,y));
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
 
ymin=@(x) 0.60*(x-1);
ymax=@(x) 1.30*(1-x);
area=integral2(nn,0,1,ymin,ymax);
 
fprintf('area = %.6f\n',area)
area = 0.958325

P3: Resolució numèrica d'EDOs

[N10] Donat el sistema d'EDOs de segon ordre

amb condicions inicials

calculeu la diferència

usant el mètode ode45. [ Ind.: passeu a un sistema de primer ordre afegint

com a funcions incògnita. ]

clear
 
P=@(t,x,y,z,w) z;
Q=@(t,x,y,z,w) w;
R=@(t,x,y,z,w) y-1.70*z+exp(t);
S=@(t,x,y,z,w) -x+2*w+1.40*sin(t);
F=@(t,X) [P(t,X(1),X(2),X(3),X(4));Q(t,X(1),X(2),X(3),X(4));...
          R(t,X(1),X(2),X(3),X(4));S(t,X(1),X(2),X(3),X(4))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[0,-1,1,2];
tf=3;
 
[t,X]=ode45(F,[t0,tf],X0,opcions);
dif=X(end,1)-X(end,2);
 
fprintf('diferencia = %.6f\n',dif)
diferencia = -326.785035

[N11] Sigui

la solució del sistema d'EDOs

amb les condicions inicials

Donada la funció

calculeu la diferència

clear
 
P=@(x,y) y-0.14*x.^2;
Q=@(x,y) -x+0.37*y^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
h=@(X) X(1).^2+X(2).^2;
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[1,0];
tf1=5;
tf2=10;
 
%usant ode45 dos cops, amb els dos temps finals tf1 i tf2:
[t,X1]=ode45(F,[t0,tf1],X0,opcions);
[t,X2]=ode45(F,[t0,tf2],X0,opcions);
dif1=h(X2(end,:))-h(X1(end,:));
fprintf('diferencia1 = %.6f\n',dif1)
diferencia1 = -0.429497
 
%usant ode45 nomes un cop, especificant tf1 i tf2:
[t,X3]=ode45(F,[t0,tf1,tf2],X0,opcions);
  %la matriu X3 nomes te 3 files, que corresponen als temps t0,tf1,tf2
dif2=h(X3(3,:))-h(X3(2,:));
fprintf('diferencia2 = %.6f\n',dif2)
diferencia2 = -0.429497

[N12] Considerem l'EDO de tercer ordre no autònoma

amb les condicions inicials

Trobeu el valor de

clear
 
R=@(t,x,y,z) 2*z-4*y-5*x+cos(t);
F=@(t,X) [X(2);X(3);R(t,X(1),X(2),X(3))];
 
opcions = odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[2.80,1.50,2.40];
tf=2;
[t,X]=ode45(F,[t0,tf],X0,opcions);
suma=X(end,1)+X(end,2)-X(end,3);
 
fprintf('suma = %.6f\n',suma)
suma = -116.386939

[N13] Donat el PVI següent per a una EDO de segon ordre no autònoma,

calculeu

aplicant el mètode ode45.

clear
 
P=@(x,y) y;
Q=@(t,x,y) -x+2*x.^2-1.30*t.*y;
F=@(t,X) [P(X(1),X(2));Q(t,X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[0.16,0.20];
tf=1;
 
[t,X]=ode45(F,[t0,tf],X0,opcions);
xf=X(end,1);
 
fprintf('x(1) = %.6f\n',xf)
x(1) = 0.274845

[N14] Donat el PVI següent per a un sistema de primer ordre no autònom,

calculeu

aplicant el mètode ode45.

clear
 
P=@(t,x,y) 0.14*t.*(x+y);
Q=@(x,y) x.^2+2*y.^2;
F=@(t,X) [P(t,X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[0.10,0.19];
tf=1;
 
[t,X]=ode45(F,[t0,tf],X0,opcions);
xf=X(end,1);
 
fprintf('x(1) = %.6f\n',xf)
x(1) = 0.126685

[N15] El moviment vertical d'un paracaigudista a l'aire ve descrit, en la hipòtesi que la fricció és proporcional al quadrat de la velocitat, per l'EDO

és l'altura respecte el terre, m és la massa, g és l'acceleració de la gravetat, i k un coeficient relacionat amb la fricció. Preneu els valors

Suposant que l'altura inicial és

i la velocitat inicial és

calculeu la distància recorreguda entre els instants

clear
 
m=83; g=9.80; k=0.30;
Q=@(x,y) -g+(k/m)*y.^2;
F=@(t,X) [X(2);Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[2830,0];
 
t1=8.70;
[t,X]=ode45(F,[t0,t1],X0,opcions);
x1=X(end,1);
 
t2=11.70;
[t,X]=ode45(F,[t0,t2],X0,opcions);
x2=X(end,1);
 
dist=x1-x2;
fprintf('distancia = %.6f\n',dist)
distancia = 149.302170

[N16] Considerem el sistema d'EDOs

i les condicions inicials

Calculeu

usant ode45.

clear
 
P=@(x,y) y.^2-0.98;
Q=@(x,y) 1.04*x.^3-1.52*y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[-1.24,-1.25];
 
tf=3;
[t,X]=ode45(F,[t0,tf],X0,opcions);
A=X(end,1);
 
tf=10;
[t,X]=ode45(F,[t0,tf],X0,opcions);
B=X(end,2);
 
suma=A+B;
fprintf('suma = %.6f\n',suma)
suma = -2.101406

[N17] Considereu el sistema d'EDOs donat pel camp

és la solució del PVI

calculeu

essent

usant ode45.

clear
 
a=-0.20; b=20.40;
 
P=@(x,y) 2*x-x.*y+a*x.*(5-x);
Q=@(x,y) -5*y+x.*y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[5,b];
 
tf=1;
[t,X]=ode45(F,[t0,tf],X0,opcions);
s1=X(end,1);
 
tf=2;
[t,X]=ode45(F,[t0,tf],X0,opcions);
s2=X(end,2);
 
solucio=s1+s2;
fprintf('solucio = %.6f\n',solucio)
solucio = 0.088523

[N18] Useu ode45 per aproximar la solució del PVI

essent

i calculeu així el valor aproximat de

clear
 
a=0.39; b=1.40;
 
P=@(x,y) sin(a*y);
Q=@(x,y) -x-y.^2/10;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[b,0];
tf=1;
 
[t,X]=ode45(F,[t0,tf],X0,opcions);
xf=X(end,1);
 
fprintf('xf = %.6f\n',xf)
xf = 1.135953

[N19] Usant ode45, resoleu l'equació del pèndol amb fricció

amb els paràmetres

i les condicions inicials

fins al temps final

Calculeu el valor de l'energia total

a l'instant final

clear
 
a=3.50; b=0.80;
Q=@(x,y) -a*sin(x)-b*y;
F=@(t,X) [X(2);Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
 
t0=0; X0=[pi/2,0];
tf=2*pi/sqrt(a);
 
[t,X]=ode45(F,[t0,tf],X0,opcions);
 
E=@(x,y) y.^2/2+a*(1-cos(x));
energia=E(X(end,1),X(end,2));
 
fprintf('energia = %.6f\n',energia)
energia = 0.219839

P4: Determinació d'esdeveniments en EDOs

[E06] Donada l'equació del pèndol amb fricció

essent

sigui

la solució que compleix les condicions inicials (angle i velocitat angular)

Trobeu l'instant

en què l'energia total

s'ha reduït a la meitat del seu valor inicial.

clear
 
g=9.80; b=0.18;
Q=@(x,y) -g*sin(x)-b*y;
F=@(t,X) [X(2);Q(X(1),X(2))];
 
E=@(x,y) y.^2/2+g*(1-cos(x));
 
t0=0; X0=[2.30,0];
E0=E(2.30,0);
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events', @(t,X) condicioE06(t,X,E,E0));
 
tf=25;
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
 
fprintf('temps = %.6f\n',te)
temps = 4.429113

[E07] Donat el sistema d'EDOs

amb condicions inicials

calculeu el valor de

al primer instant

en què la solució torna a tallar la semirecta

clear
 
P=@(x,y) 3*x-2*y+0.30*sin(x.*y);
Q=@(x,y) 4*x-y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE07);
 
t0=0; X0=[0.80,0.80];
tf=5;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
jtall=find( Xe(:,1)>0 );
xtall=Xe(jtall(2),1); %notem que jtall(1) correspon a la condicio inicial
 
fprintf('xtall = %.6f\n',xtall)
xtall = 22.151235

[E08] Considereu l'EDO

essent

que correspon al moviment d'un pèndol sense fricció de massa

Sabem que l'energia

és una quantitat conservada. Calculeu el període de l'òrbita periòdica corresponent al nivell d'energia

[ Ind.: heu d'escollir una condició inicial sobre aquesta corba de nivell, per exemple trobant la seva intersecció amb

]

clear
 
a=14.50;
Q=@(x,y) -a*sin(x);
F=@(t,X) [X(2);Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE08);
 
E0=1.70;
y0=sqrt(2*E0); %aillant y=+sqrt(...) de E(0,x')=E0 i avaluant en x=0
t0=0; X0=[0,y0];
tf=25;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
temps=te(3); %volem el segon tall
 
fprintf('periode = %.6f\n',temps)
periode = 1.675058

[E09] Considerem l'EDO lineal de segon ordre

que correspon a un moviment oscil·latori amb fricció donada pel coeficient

i sobre el qual actua una força externa donada per la funció

Sigui

i suposem que

(per tant la força externa només actua a partir de l'instant

). Suposant condicions inicials

trobeu el primer instant

en què s'assoleix la posició

clear
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events', @condicioE09);
 
%primera opcio:
%  definim una primera EDO amb f(t)=0 i la resolem fins a l'instant T0,
%  i prenem el punt final X1(end,:) obtingut, com a condicio inicial
%  d'una segona EDO amb f(t)=t-T0
 
Q1=@(x,y) -8*x-0.19*y;
F1=@(t,X) [X(2);Q1(X(1),X(2))];
 
T0=10;
[t1,X1]=ode45(F1,[0,T0],[0,0.70],opcions);
 
f=@(t) t-T0;
Q2=@(t,x,y) -8*x-0.19*y+f(t);
F2=@(t,X) [X(2);Q2(t,X(1),X(2))];
 
tf=T0+20;
[t2,X2,te,Xe]=ode45(F2,[T0,tf],X1(end,:),opcions);
fprintf('te = %.6f\n',te)
te = 17.883356
 
%segona opcio:
%  definim una unica EDO amb una funcio fb(t) que inclou els dos casos
 
fb=@(t) max(0,t-T0);
Qb=@(t,x,y) -8*x-0.19*y+fb(t);
Fb=@(t,X) [X(2);Qb(t,X(1),X(2))];
 
[tb,Xb,teb,Xeb]=ode45(Fb,[0,30],[0,0.70],opcions);
fprintf('teb = %.6f\n',teb)
teb = 17.883356

[E10] La solució del sistema

amb condicions inicials

talla la circumferència

en almenys dos punts. Trobeu l'instant

en què té lloc el segon tall.

clear
 
P=@(x,y) y-2.13*x.*y;
Q=@(x,y) -x+sin(y);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events', @condicioE10);
 
t0=0; X0=[1.59,1.50];
tf=3;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
fprintf('t2 = %.6f\n',te(2))
t2 = 1.542344

[E11] La solució del PVI

és periòdica. Calculeu-ne el període.

clear
 
P=@(x,y) -x.*y;
Q=@(x,y) (x-1.20).*(y+1.90);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
t0=0; X0=[2.10,0];
DX0=F(0,X0);
Dy0=DX0(2); %com que es positiu, posarem signe=1 a la funcio condicio
fprintf('Dy0 = %.6f\n',Dy0)
Dy0 = 1.710000
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE11);
 
tf=25;
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
 
periode=te(2); %si fem signe=0, el periode sera te(3)
fprintf('periode = %.6f\n',periode)
periode = 4.269537

[E12] Sigui

la solució del sistema

amb condicions inicials

és l'instant en què la coordenada

assoleix el seu primer màxim, trobeu el valor de

[ Ind.: en un màxim, la derivada

passa de positiva a negativa. ]

clear
 
P=@(x,y) y+x.^2+0.11*x.^3-y.^2;
Q=@(x,y) -x-2*x.*y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events', @(t,X) condicioE12(t,X,Q));
 
t0=0; X0=[0,-0.34];
tf=5;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
y1=Xe(1,2);
 
fprintf('y1 = %.6f\n',y1)
y1 = 0.462813

[E13] Donada l'equació del pèndol amb fricció

essent

sigui

la solució que compleix les condicions inicials (angle i velocitat angular)

(l'angle inicial correspon a la posició vertical superior). La velocitat inicial és suficient perquè el pèndol doni més d'una volta completa. Quina és la velocitat angular (mesurada a la posició vertical superior) quan ha donat 2 voltes completes? [ Ind.: l'angle s'incrementa en

després de cada volta completa. ]

clear
 
g=9.80; b=0.10;
P=@(x,y) y;
Q=@(x,y) -g*sin(x)-b*y;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE13);
 
t0=0; X0=[-pi,6.70];
tf=10;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
vf=Xe(1,2);
 
fprintf('vf = %.6f\n',vf)
vf = 5.141403

[E14] Sigui el sistema d'EDOs

i preneu la condició inicial

Trobeu el temps necessari per assolir la recta

amb

clear
 
P=@(x,y) y;
Q=@(x,y) -x+x.^2;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE14);
 
t0=0; X0=[1.02,-0.13];
tf=7;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
temps=te(1);
 
fprintf('temps = %.6f\n',temps)
temps = 3.997270

[E15] Sigui

la solució del sistema

amb les condicions inicials

Trobeu el temps que triga en tornar a tallar el semieix positiu de les

clear
 
P=@(x,y) 3.50*(x+4*y)+x.^2+4*x.*y;
Q=@(x,y) -6*x-5*y+(6*x.*y+5*y.^2)/3.30;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE15);
 
t0=0; X0=[2.30,0];
tf=20;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
j=1;
abscissa=-1;
while( abscissa<=0 )
  j=j+1;
  temps=te(j);
  abscissa=Xe(j,1);
end
 
fprintf('temps = %.6f\n',temps)
temps = 0.768799

[E16] La solució de l'EDO de segon ordre

amb condicions inicials

defineix una funció

que assoleix un màxim relatiu en un instant

Trobeu el valor de

clear
 
P=@(x,y) y;
Q=@(x,y) -7.37*y-5.81-6.02*sin(x);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE16);
 
t0=0; X0=[1.14,2.98];
tf=1;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
xmax=Xe(1);
 
fprintf('xmax = %.6f\n',xmax)
xmax = 1.318213

[E17] Sigui

la solució del sistema lineal homogeni atractor

essent

amb les condicions inicials

Trobeu el primer instant t en què la funció

assoleix la meitat del seu valor inicial.

clear
 
a=1.10; b=1.80;
A = [-1,a,b; 0,-1,1; 0,0,-3];
F=@(t,X) A*X;
  % Nota: Tambe podem definir les 3 components per separat,
  %   P=@(x,y,z) -x+a*y+b*z;
  %   Q=@(x,y,z) -y+z;
  %   R=@(x,y,z) -3*z;
  %   F=@(t,X) [P(X(1),X(2),X(3));Q(X(1),X(2),X(3));R(X(1),X(2),X(3))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE17);
 
t0=0; X0=[1,1,1];
tf=4;
 
[t,X,te,Xe] = ode45(F, [t0,tf], X0, opcions);
temps=te(1);
 
fprintf('temps = %.6f\n',temps)
temps = 2.417505

[E18] Considerem el pèndol sense fricció modelitzat com el sistema d'EDOs

Considerem, com a condició inicial, el punt

que pertany a la solució asimptòtica entre els punts d'equilibri

i continguda al semiplà

(heu d'escollir

de manera que el valor de la quantitat conservada

en aquest punt coincideixi amb el valor a

). Quant de temps triga la solució a tallar la recta

? Preneu els valors

clear
 
g=9.80; l=0.54; %no necessitem m
Q=@(x,y) -(g/l)*sin(x);
F=@(t,X) [X(2);Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE18);
 
ymax=2*sqrt(g/l);
  %obtingut aillant y=+sqrt(...) de la quantitat conservada,
  %i avaluant en x=0
t0=0; X0=[0,ymax];
tf=25;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
 
temps=te(1);
fprintf('temps = %.6f\n',temps)
temps = 0.094646

[E19] Sigui

la solució del sistema d'EDOs

amb condicions inicials

essent

és el primer instant en què la solució talla la paràbola

essent

trobeu el valor de

clear
 
a=0.30;
 
P=@(x,y) -sin(x)+y;
Q=@(x,y) cos(x.*y);
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
 
opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
opcions=odeset(opcions,'Events',@condicioE19);
 
t0=0; X0=[0,a];
tf=5;
 
[t,X,te,Xe]=ode45(F,[t0,tf],X0,opcions);
yt1=Xe(1,2);
 
fprintf('yt1 = %.6f\n',yt1)
yt1 = 1.506721

Funcions auxiliars per a P4

function [expr,aturar,signe]=condicioE06(t,X,E,E0)
  expr=E(X(1),X(2))-E0/2;
  aturar=0;
  signe=-1;
end
 
function [expr,aturar,signe]=condicioE07(t,X)
  expr=X(2)-X(1);
  aturar=0;
  signe=0;
end
 
function [expr,aturar,signe]=condicioE08(t,X)
  expr=X(1);
  aturar=0;
  signe=0;
end
 
function [expr,aturar,signe]=condicioE09(t,X)
  expr=X(1)-1;
  aturar=1;
  signe=0;
end
 
function [expr,aturar,signe]=condicioE10(t,X)
  expr=X(1)^2+X(2)^2-4;
  aturar=0;
  signe=0;
end
 
function [expr,aturar,signe]=condicioE11(t,X)
  expr=X(2);
  aturar=0;
  signe=1;
end
 
function [expr,aturar,signe]=condicioE12(t,X,Q)
  expr=Q(X(1),X(2));
  aturar=0;
  signe=-1;
end
 
function [expr,aturar,signe]=condicioE13(t,X)
  expr=X(1)-3*pi;
  aturar=0;
  signe=1;
end
 
function [expr,aturar,signe]=condicioE14(t,X)
  expr=X(2)+0.41;
  aturar=1;
  signe=1;
end
 
function [expr,aturar,signe]=condicioE15(t,X)
  expr=X(2);
  aturar=0;
  signe=0;
end
 
function [expr,aturar,signe]=condicioE16(t,X)
  expr=X(2);
  aturar=0;
  signe=-1;
end
 
function [expr,aturar,signe]=condicioE17(t,X)
  expr=X(1)-1/2;
  aturar=1;
  signe=-1;
end
 
function [expr,aturar,signe]=condicioE18(t,X)
  expr=X(1)-pi/4;
  aturar=1;
  signe=1;
end
 
function [expr,aturar,signe]=condicioE19(t,X)
  b=1.20;
  expr=X(2)-b*X(1)^2;
  aturar=0;
  signe=0;
end