Pràctica: RESOLUCIÓ NUMÈRICA D'EDOs

Veurem com podem integrar (és a dir, resoldre) numèricament un sistema d'equacions diferencials ordinàries (EDOs). Començarem pel cas d'una sola EDO, amb el mètode d'Euler i el mètode RK45F. Després ho veurem per a un sistema d'EDOs.

Mètode d'Euler

Plantegem un problema de valor inicial (PVI) per a una EDO de primer ordre,

on f és una funció escalar de 2 variables i

són valors donats. Volem trobar aproximacions de la solució

per a valors de t en un interval

El mètode més simple per resoldre numèricament una EDO és el mètode d'Euler, que es basa en la idea següent. Escollim un pas d'integració

i considerem valors equidistants

Per obtenir una aproximació

usem:

i per tant és raonable considerar

D'aquesta manera podem anar calculant cada aproximació

a partir de l'anterior, aplicant les iteracions següents:

Es diu que el mètode d'Euler és d' ordre 1 , ja que es basa en una aproximació de la derivada que té un error

De fet, es pot veure que l'aproximació del valor final

també té un error

(com es pot comprovar a l'exemple de més avall).

Al fitxer odeEuler.m podeu trobar una funció per aplicar el mètode d'Euler a un PVI amb un pas d'integració donat: odeEuler(f, [t0,tf], x0, h). (Nota: També podem integrar "enrere" si escollim

)

Com a exemple, volem aproximar la solució del PVI

per a

amb passos d'integració

Dibuixarem les gràfiques de la solució aproximada obtinguda en cada cas, i escriurem el valor aproximat de

També compararem amb la solució exacta

close all %tanquem figures anteriors (optatiu)

figure %nova figura

hold on

f=@(t,x) (t-1)/2.*x; %EDO a resoldre

t0=0; x0=3; %condicions inicials

tf=1; %temps final

for h=[0.2,0.1,0.05] %passos d'integracio

[t,x]=odeEuler(f,[t0,tf],x0,h);

%retorna els valors de t i x com dos vectors columna

plot(t,x,'b') %grafica formada per segments

plot(t,x,'ro') %marquem els punts entre cada dos segments

fprintf('h=%.2f, x(%.2f) ~ %f\n',h,tf,x(end))

%valor aproximat de x(tf) amb aquest pas h

%(es el darrer element del vector columna x)

end

h=0.20, x(1.00) ~ 2.196730

h=0.10, x(1.00) ~ 2.267481

h=0.05, x(1.00) ~ 2.302144

xlabel('t') %nom abscissa

ylabel('x') %nom ordenada

xx=@(t) 3*exp((t.^2-2*t)/4); %solucio exacta

t=0:0.01:1;

plot(t,xx(t),'k')

fprintf('x(%.2f) = %f\n',tf,xx(tf)) %valor exacte de x(tf)

x(1.00) = 2.336402

Mètode RK45F: la funció ode45

Per millorar les aproximacions, cal usar mètodes d'ordre més alt que el d'Euler, a base d'avaluar f en diversos punts adequadament escollits, en cada iteració. Tenim la funció de Matlab ode45 , corresponent a una variant del mètode de Runge-Kutta-Fehlberg d'ordres 4 i 5 ( RK45F ), que és una combinació de mètodes d'ordres 4 i 5 que permet un control automàtic del pas d'integració (que anirà variant en cada iteració) a partir d'una tolerància d'error donada. Així doncs, no hem de donar cap pas d'integració sinó fixar prèviament la magnitud de l'error que volem permetre.

Il·lustrem-ho amb l'EDO considerada anteriorment. En aquest cas, dibuixarem la solució obtinguda per a

amb 5 condicions inicials diferents:

figure %nova figura

hold on

opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);

f=@(t,x) (t-1)/2.*x; %EDO a resoldre

t0=0; %temps inicial

tf=3; %temps final

for x0=2:0.5:4 %valors inicials

[t,x]=ode45(f,[t0,tf],x0,opcions);

plot(t,x)

end

xlabel('t')

ylabel('x')

Nota: Quan es resol un PVI per a una EDO no lineal, segons l'interval d'integració escollit pot passar que no es pugui trobar la solució a tot l'interval perquè s'apropi a infinit per algun

(és a dir, que la gràfica de

tingui una asímptota vertical). En aquest cas, ens apareixeran missatges del tipus "Warning: Failure at t=(...). Unable to meet integration tolerances without reducing the step size below the smallest value allowed (...) at time t". Per exemple, ho podem veure si intentem resoldre el PVI

a l'interval

La solució

només té sentit per a

Resolució d'un sistema d'EDOs amb ode45

De manera semblant, podem usar la funció ode45 per aproximar la solució d'un PVI per a un sistema d'EDOs. Com a exemple, per al sistema 2D (autònom)

amb les condicions inicials

trobarem la solució aproximada per a

i veurem diferents maneres de representar-la gràficament.

opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);
P=@(x,y) y-x.^3; %components del camp vectorial
Q=@(x,y) -x;
F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];
  % Nota 1:
  %   definim F=[P;Q] com un vector COLUMNA, que ha de dependre
  %   de la variable independent t (encara que sigui autonom)
  %   i de les coordenades X=[x,y]
  % Nota 2:
  %   tambe podem definir F directament, sense definir previament
  %   les seves components P i Q, escrivint
  %   "F=@(t,X) [X(2)-X(1).^3; -X(1)];"
t0=0; X0=[1,0]; %condicions inicials
tf=20; %temps final
[t,X]=ode45(F,[t0,tf],X0,opcions);
  %retorna els valors de t com un vector columna,
  %i els valors de X(t)=(x(t),y(t)) com una matriu de 2 columnes

Per representar gràficament la solució tenim diferents opcions. Representarem

com dues gràfiques, i després

com una corba parametritzada.

Ho inclourem tot en una "matriu de gràfiques" 2x2, amb la columna de l'esquerra (elements 1 i 3) ocupada per les gràfiques

i tota la columna de la dreta (elements [2,4]) per la representació de

figure %nova figura

hold on

subplot(2,2,1)

plot(t,X(:,1)) %grafica de x(t)

xlabel('t'), ylabel('x')

subplot(2,2,3)

plot(t,X(:,2)) %grafica de y(t)

xlabel('t'), ylabel('y')

subplot(2,2,[2,4])

plot(X(:,1),X(:,2)); %corba (x(t),y(t))

axis equal

xlabel('x'), ylabel('y')

Podem calcular de forma successiva moltes solucions amb condicions inicials diferents. En el mateix exemple anterior, calcularem la solució amb condicions inicials

per a

figure %nova figura

hold on

opcions=odeset('AbsTol',1e-8,'RelTol',1e-8);

P=@(x,y) y-x.^3;

Q=@(x,y) -x;

F=@(t,X) [P(X(1),X(2));Q(X(1),X(2))];

t0=0; tf=20;

X0=[1,0; 2,0; 3,0]; %matriu de condicions inicials (una a cada fila)

nci=size(X0,1); %nombre de condicions inicials

for i=1:nci

[t,X]=ode45(F,[t0,tf],X0(i,:),opcions);

plot(X(:,1),X(:,2))

end

axis equal

xlabel('x'), ylabel('y')

________