EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
PONTOS DE EQUILÍBRIO
Consideremos o sistema autônomo
dy
= f ( y ) . O ponto de equilíbrio é definido pela soluções de
dt
f ( y) = 0 .
Exemplo 1:
dy
= y (1 − y )
dt
O arquivo .m derivy
function dy = derivy(t,y)
dy=y.*(1-y);
é chamado do workspace do MATLAB:
%ponto de equilibrio
y=linspace(-0.5,1.5,100);
plot(y,derivy(t,y))
para produzir:
As soluções de equilíbrio cortam o plano y em REGIÕES INDEPENDENTES.
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
O gráfico acima foi produzido pelo seguinte arquivo .m
%ponto de equilibrio
y=linspace(-0.5,1.5,100);t=[];
figure(1)
plot(y,derivy(t,y))
%ponto inicial x0
x0=linspace(-0.5,1.5,6);
cor=['b';'m';'g';'y';'r';'k'];
tspan=[0 5];
figure(2), hold on
for i=1:length(x0)
[t,y]=ode45('derivy',tspan',x0(i));
plot(t,y,cor(i)),plot(t(1),y(1),[cor(i) 'o']);
pause
end
[t,y]=ode45('derivy',tspan',1);
plot(t,y,cor(i)),plot(t(1),y(1),'k--');
[t,y]=ode45('derivy',tspan',0);
plot(t,y,cor(i)),plot(t(1),y(1),'k--');
axis([0 1 -2 2])
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Note-se que:
1)
y(0)<0 ⇒ f(y)<0 ⇒ y sempre cai
0<y(0)<1 ⇒ f(y)>0 ⇒ y sempre aumenta
y(0)>1 ⇒ f(y)<0 ⇒ y sempre cai
2)
O ponto de equilíbrio y=1 “atrai” trajetórias enquanto o ponto y=0 “repele” trajetórias.
Exemplo 2
dy
= 0.2 y (5 − y )( y − 2)
dt
function dy = derivy2(t,y)
dy=0.2*y.*(5-y).*(y-2);
%ponto de equilibrio
y=linspace(-1,7,100);t=[];
figure(1)
plot(y,derivy2(t,y))
%ponto inicial x0
x0=linspace(-1,7,6);
cor=['b';'m';'g';'y';'r';'k'];
tspan=[0 5];
figure(2), hold on
for i=1:length(x0)
[t,y]=ode45('derivy2',tspan',x0(i));
plot(t,y,cor(i)),plot(t(1),y(1),[cor(i) 'o']);
pause
end
axis([0 1 -2 7])
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
A classificação dos pontos de equilíbrio está apresentada no gráfico a seguir.
Convenção utilizada:
f(y)>0 ↑
f(y)<0 ↓
Exemplo 3
dy1
= y1 (1 − y1 ) − y1 y 2
dt
dy 2
= 2 y 2 (1 − y 2 / 2) − 3 y1 y 2
dt
Pontos de equilíbrio
dy1
= 0 ⇒ y1 (1 − y1 ) − y1 y 2 = 0 ⇒ 1 − y1 − y 2 = 0 ou y1 = 0 (isóclina 1)
dt
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
dy 2
= 0 ⇒ 2 y 2 (1 − y 2 / 2) − 3 y1 y 2 = 0 ⇒ 2 − 3 y1 − y 2 = 0 ou y 2 = 0 (isóclina 2)
dt
Construindo o PLANO DE FASE (y1 x y2) com o MATLAB:
%ponto de equilibrio
tspan=[0 5];
hold on
y10=[0.1 0.3 0.4 0.6 1 1.5];
y20=[0 0.8 1.0 1.5 2 2.5];
cor=['b';'m';'g';'y';'r';'k'];
for i=1:length(y10)
for j=1:length(y20)
[t,y]=ode45('derivy3',tspan',[y10(i) y20(j)]);
plot(y(:,1),y(:,2),cor(i)),plot(y(1,1),y(1,2),[cor(i) 'o']);
pause
end
end
axis([0 1.5 0 2.5])
%isoclinas
y1=linspace(0,1.2,12)';
isoc1=[y1 1-y1];
isoc2=[y1 2-3*y1];
isoc3=[zeros(size(y1)) y1];
isoc4=[y1 zeros(size(y1))];
plot(isoc1(:,1),isoc1(:,2),'c')
plot(isoc2(:,1),isoc2(:,2),'c')
plot(isoc3(:,1),isoc3(:,2),'c')
plot(isoc4(:,1),isoc4(:,2),'c')
function dy = derivy3(t,y)
dy(1)=y(1)*(1-y(1))-y(1)*y(2);
dy(2)=2*y(2)*(1-y(2)/2)-3*y(1)*y(2);
dy=dy(:);
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
PLANO DE FASES
Dado o sistema autônomo
dx
= Ax
dt
Ax = λ x
assumindo que A tem uma base de n autovetores x1...xn com autovalores λ1... λn
y = c1e λ1t + ... + c n e λnt
λ1... λn devem ter parte real positiva para que o sistema seja estável.
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Exemplo 4
2 1 
dx
= Ax , A = 

dt
2 − 1
>> [V,D]=eig([2 1;2 -1])
V =
0.8719
0.4896
-0.2703
0.9628
2.5616
0
0
-1.5616
D =
onde
λ
D= 1
0
0
λ 2 
λ1>0 ⇒ V1 é um subespaço instável
λ2<0 ⇒ V2 é um subespaço estável
Ponto de equilíbrio: A x = 0
>> A=[2 1;2 -1]
A =
2
2
>> b=[0;0]
b =
0
0
1
-1
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
>> x=A\b
x =
0
0
%ponto de equilibrio
tspan=[0 10];
hold on
y10=[-2 -1 0 1 2];
y20=[-2 -1 0 1 2];
cor=['b';'m';'g';'y';'r';'k'];
for i=1:length(y10)
for j=1:length(y20)
[t,y]=ode45('derivy4',tspan',[y10(i) y20(j)]);
plot(y(:,1),y(:,2),cor(i)),plot(y(1,1),y(1,2),[cor(i) 'o']);
pause
end
end
axis([-2 2 -2 2])
function dy = derivy4(t,x)
A=[2 1;2 -1];
dy=A*x;
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Exemplo 4b
− 3 1 
A=

 1 − 3
>> [V,D]=eig([-3 1;1 -3])
V =
0.7071
-0.7071
0.7071
0.7071
D =
-4
0
0
-2
Sistema estável (2 autovalores negativos )
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
function dy = derivy4b(t,x)
A=[-3 1;1 -3];
dy=A*x;
%ponto de equilibrio
tspan=[0 10];
hold on
y10=[-2 -1 0 1 2];
y20=[-2 -1 0 1 2];
cor=['b';'m';'g';'y';'r';'k'];
for i=1:length(y10)
for j=1:length(y20)
[t,y]=ode45('derivy4b',tspan',[y10(i) y20(j)]);
plot(y(:,1),y(:,2),cor(i)),plot(y(1,1),y(1,2),[cor(i) 'o']);
pause
end
end
axis([-2 2 -2 2])
Neste exemplo, ocorre um nó estável ou “próprio” (“sink”)
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Exemplo 4c:
1 0
A=

0 1 
>> [V,D]=eig(eye(2,2))
V =
1
0
0
1
1
0
0
1
D =
Neste exemplo, ocorre um nó instável ou “impróprio” (“source”)
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Exemplo 4d:
1 0 
A=

1 − 1
>>
[V,D]=eig([1 0; 1 -1])
V =
0
1.0000
0.8944
0.4472
D =
-1
0
0
1
Neste caso, ocorre um “ponto de sela”.
EQE-709 – SIMULAÇÃO E CONTROLE DE PROCESSOS
Profa. Ofélia de Q.F. Araújo - Escola de Quimica / UFRJ
Exemplo 4e:
 0 1
A=

 − 4 0
>> [V,D]=eig([0 1; -4 0])
V =
0 - 0.4472i
0.8944
0 + 0.4472i
0.8944
0 + 2.0000i
0
0
0 - 2.0000i
D =
Para autovalores complexos, puramente imaginários, ocorrem trajetórias fechadas, e o ponto de
equilíbrio é CENTRAL.