MATLAB Básico
Carlos André Vaz Junior
[email protected]
http://www.eq.ufrj.br/links/h2cin/carlosandre
MATLAB Básico
Mais de 1 milhão de resultados
MATLAB Básico
?
http://newsreader.mathworks.com
MATLAB Básico
MATLAB Básico
MATLAB Básico
MATLAB Básico
MATLAB Básico
Agora a = 2, faço tudo de novo?!
MATLAB Básico
MATLAB Básico
MATLAB Básico
Tipos Básicos
MATLAB Básico
Matriz
Case Sensitive!
Char Array
Estrutura
CaSe SeNsItIvE!
MATLAB Básico
Criando uma matriz:
MATLAB Básico
Criando um “char array”:
Banco de Dados da “Turma”:
Alunos: Carla, João, Bruno, Luis, Marcela
Professor: Marcelo
Horário: 13h
Sala: 221
MATLAB Básico
Estrutura:
turma.alunos.nomes=strvcat( 'carla',’joao','bruno', ...
'luis', 'marcela‘ );
turma.professor.nome=(‘Marcelo‘)
turma.horario=1300
turma.sala=221
MATLAB Básico
Comando “who” e “whos”
Use “;” para evitar que o resultado apareça na tela.
Use A=0:0.5:10 para gerar matrizes com dados em seqüência.
Use “clear A” para apagar a variável A.
Use “clear all” para apagar todas as variáveis armazenadas.
MATLAB Básico
Use “size(A) ” para identificar as dimensões da matriz. A maior
dimensão é dada pelo comando “length(A) ”
Dicas!
i) Soma e subtração: soma (ou subtrai) elemento por elemento da
matriz.
A+B
A-B
MATLAB Básico
ii) Multiplicação e Divisão de matrizes: atenção às regras da álgebra,
pois as dimensões das matrizes têm que ser coerentes!
A*B
A/B
iii) Multiplicação e divisão elemento por elemento:
A .* B
A ./ B
iv) Matriz Transposta:
A’
v) Cria Matriz Identidade:
eye(número de linhas, número de colunas)
vi) Cria Matriz Zeros:
zeros(número de linhas, número de colunas)
MATLAB Básico
vii) Cria Matriz Uns:
ones(número de linhas, número de colunas)
viii) Cria Matriz Randômica (composta de números aleatórios):
rand(número de linhas, número de colunas)
ix) Determinante:
det(matriz)
x) Inversa:
inv(matriz)
MATLAB Básico
xi) Dimensões da matriz:
size(matriz)
lenght(matriz)
numel(matriz)
Veja também: flipud e fliplr
MATLAB Básico
1
5
9
13
2
6
10
14
3
7
11
15
4
8
12
16
Elemento = Matriz(2,3) ou Matriz(10)
MATLAB Básico
Programa Principal / Workspace
global C
C=100
D=22
MATLAB Básico
Função Alfa
A=1
B=2
global C
C=100
Função Beta
E=15
F=55
C=23
Achando a posição do menor valor de uma matriz:
x=[1 2 3 4 5 6;
2 1 3 3 2 1];
MATLAB Básico
%Forma linear:
xmin=min(x);
xmin=min(xmin);
[i,j]=find(x==xmin);
%Forma condensada:
[i,j]=find(x==(min(min(x))));
Achando o menor valor de uma função:
X = fzero('sin',2)
MATLAB Básico
função
estimativa inicial
Veja também: fsolve e fmin
MATLAB Básico
if:
MATLAB Básico
AND
OR
Falso
Verdadeiro
MATLAB Básico
AND
a
b
1
1
0
OR
resultado
a
b
resultado
1
1
1
1
1
0
0
1
1
1
0
0
1
0
1
0
0
0
0
0
0
MATLAB Básico
Case:
switch I
case 1,
disp('I vale 1')
case 2,
disp('I vale 2')
otherwise
disp('I nao eh nem 1 nem 2')
end
While:
MATLAB Básico
while I < 10,
disp(‘oi’);
I=I+1;
end
Manipule o ponteiro I na rotina executada
pelo “while”
For:
MATLAB Básico
for J = 1:100,
A(1,J) = 1/(I+J-1);
end
Incremento automático do ponteiro J
a cada loop
MATLAB Básico
>> figure(1)
>> figure(2)
>> t=0:0.01:10;
>> y=sin(t);
>> plot(t,y)
>> z=cos(t);
>> plot(t,z)
Use “[x,y]=ginput(2)” para capturar dois pontos no gráfico
Use “close all” para fechar todas as figuras
MATLAB Básico
Use “clf” para apagar a figura atual
Dica!
MATLAB Básico
>> figure(3)
>> plot(t,y)
>> subplot(1,2,1)
>> subplot(1,2,2)
>> plot(t,z)
MATLAB Básico
>>
>>
>>
>>
>>
>>
>>
t=0:0.25:10;
y=sin(t);
plot(t,y,'r+')
xlabel('tempo')
ylabel('seno')
title('Seno vs. Tempo')
Axis([0 10 -2 2])
>>
>>
>>
>>
>>
t=0:0.01:10;
y=sin(t);
z=cos(t);
plot(t,y,'g-',t,z,'r-')
legend('seno','cosseno')
MATLAB Básico
Ou...
>> t=0:0.01:10;
>> y=sin(t);
>> z=cos(t);
>> plot(t,y,'g-‘)
>> hold on
>> plot(t,z,'r-')
>> legend('seno','cosseno')
MATLAB Básico
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
Z=exp(-0.5*(X.^2+Y.^2));
MATLAB Básico
colormap jet
figure(1);surf(X,Y,Z); rotate3d on; shading interp;
MATLAB Básico
MATLAB Básico
MATLAB Básico
%Malha triangular da base
%malha da base
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
Z=1-X-Y;
%aplica a restrição para usar só a base do triangulo
%onde existe consistência física (o que nao tem vira "Not a Number")
iz=find(Z<0);Z(iz)=nan;
Composição
(3 componentes)
MATLAB Básico
%Malha triangular da base
%malha da base
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
Z=1-X-Y;
%aplica a restrição para usar só a base do triangulo
%onde existe consistência física (o que não tem vira "Not a Number")
iz=find(Z<0);Z(iz)=nan;
Alguns Z são
negativos! Não pode!
MATLAB Básico
%Malha triangular da base
%malha da base
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
Z=1-X-Y;
%aplica a restrição para usar só a base do triangulo
%onde existe consistência física (o que não tem vira "Not a Number")
iz=find(Z<0);Z(iz)=nan;
vv1=(X.*log(X))+(Y.*log(Y))+(Z.*log(Z));
MATLAB Básico
%gráfico da superfície
colormap jet
figure(1);surf(X,Y,vv1); rotate3d on; shading interp;
xlabel('X1');ylabel('X2');zlabel('DeltaGi/RT');
MATLAB Básico
MATLAB Básico
Exemplos
MATLAB Básico
Exemplo
1
MATLAB Básico
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
pode-se escrever o balanço de massa do sistema
Ainda,
dmt 
  FE  F 
dt
(2)
dmt 
dht 
 A
dt
dt
(3)
dht  1
 FE  F 
dt
A
(4)
MATLAB Básico
e, portanto,
1
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
Freqüentemente, considera-se a vazão de saída do tanque
proporcional
à
altura
da
coluna
de
líquido
é inversamente
proporcional a uma resistência ao escoamento (R):
F 
h
R
(5)
MATLAB Básico
Logo,
dht  1 
h
  FE  
dt
A
R
(6)
1
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
Este modelo simples de um tanque de nível, sem balanço de
energia, possui uma solução analítica:
t



ht   RFE 1  e RA 


(7)
Para simular este modelo, basta escolher os valores das
MATLAB Básico
constantes R, A e FE, das condições iniciais h0 e t0.
A simulação da solução analítica do modelo do tanque de
nível é mostrada a seguir.
1
MATLAB Básico
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10;
% m3/h
% Tempo de simulação
t = 0.0 : 0.01 : 10.0;
%h
% Simulação da altura de líquido
h = R*Fe*(1 - exp(-t/(R*A)));
%m
% Visualização da simulação
plot(t,h);
title('Simulação do tanque de nível');
xlabel('Tempo (h)');
ylabel('Altura (m)');
MATLAB Básico
Verifique a consistência do calculo: a
matriz “h” gerada também deve ser 1x1000, já que
cada instante “t” gerou um valor “h”. É sempre útil
conferir a dimensão das variáveis, principalmente a
medida que as rotinas forem tornando-se
complexas.
Dica!
MATLAB Básico
Exemplo
2
Modelagem & Dinâmica de Processos
Muitas vezes é muito trabalhoso, ou mesmo impossível,
Modelos simples - o tanque de nível
encontrar a solução analítica para o conjunto de
equaçõesFreqüentemente,
diferenciais. Nesse
casoatemos
quesaída
simular
considera-se
vazão de
do tanque
usando
soluçãoà numérica
equações
diferenciais.
proporcional
altura da das
coluna
de líquido
é inversamente
Vamos
assumir aque
modelo ao
doescoamento
exemplo 1
não tivesse
proporcional
umaoresistência
(R):
solução analítica, e então usar o Matlab para estudar o
comportamento da altura dohnível com o tempo. A
F 
(5)
equação diferencial
será:
R
MATLAB Básico
Logo,
dht  1 
h
  FE  
dt
A
R
(6)
1
MATLAB Básico
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10; % m3/h
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; % h
% Simulação da altura de líquido
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
% Visualização da simulação
plot(t,h);
title('Simulação do tanque de nível');
xlabel('Tempo (h)');
ylabel('Altura (m)');
function dh = dhdt(t,h,flag,par)
R = par(1);
A = par(2);
Fe = par(3);
dh = (Fe-(h/R))/A;
Nesse caso temos uma equação diferencial, então deveremos
usar uma função Matlab específica para a resolução de eq.
diferenciais. No caso temos a ODE45. A função ODE45
implementa um esquema de solução de sistemas de EDO’s por
método de Runge-Kutta de ordem média (consulte o help sobre
ODE45 para maiores detalhes).
MATLAB Básico
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
Os parâmetros enviados entre parênteses são aqueles que
devemos passar para a ODE45:
-1º argumento de ode45 é uma string contendo o nome do arquivo .m
com as equações diferenciais. Neste caso, o arquivo chama-se dhdt.m.
MATLAB Básico
-2º argumento é um vetor que pode conter (i) dois elementos: os
tempos inicial e final da integração, ou (ii) todos os valores de tempo para
os quais deseja-se conhecer o valor da variável integrada.
-3º argumento é o vetor contendo as condições iniciais das variáveis
dependentes das EDO’s. Os valores dos elementos do vetor de condições
iniciais precisam estar na mesma ordem em que as variáveis
correspondentes são calculadas na função passada como 1º argumento
para ode45 (neste caso, dhdt.m). Nesse caso em particular só temos uma
variável dependente, assim temos uma única condição inicial.
-4º argumento é o vetor de opções de ode45. Há várias opções do
método que podem ser ajustadas. Entretanto, não deseja-se alterar
os valores-padrão. Neste caso, é passado um vetor vazio, apenas para
marcar o lugar das opções.
MATLAB Básico
-5º argumento é um vetor contendo parâmetros de entrada para a
função dhdt.m. Observe que a função .m deve ler esses parâmetros
na ordem correta (recebe como variável local “par”).
Os resultados da simulação são obtidos nos dois parâmetros
entre colchetes (t , h).
A codificação do arquivo .m segue o mesmo formato já explicado para
funções porém com algumas particularidades.
No caso específico de um arquivo .m que deve ser chamado por uma
função de solução EDO’s (todas as ODExx), a declaração deste arquivo
deve seguir a sintaxe:
MATLAB Básico
function dy = nomefun(t, y, flag, arg1, ..., argN)
onde
•dy é o valor da(s) derivada(s) retornadas
•t e y são as variáveis independente e dependente, respectivamente.
•Opcional: caso deseje-se receber outros parâmetros, a função deve
receber um argumento marcador de lugar chamado flag. Após este,
ela recebe quaisquer outros parâmetros.
MATLAB Básico
Exemplo
3
MATLAB Básico
Modelagem & Dinâmica de Processos
Modelagem & Dinâmica de Processos
Modelos simples - tanque de aquecimento
Modelos simples - tanque de aquecimento
Como no caso anterior, o balanço de massa pode ser escrito
   V dT
d :VT
como
dt
dt
T
dV
dh 
 dT
 A h
T

dt
dt 
 h dt
1
dht 
  FE  
dt
A
R
 dT T 
h 
C p A  h
  FE     FE HE  FH   Q
dt
A
R
MATLAB Básico

 como:
O
balanço de energia
é escrito
 UTFh H FH
dT C1 d FVT
FE   Q
U 
E TE  
 p


E 
E T
 A

dt
h  dtA
C p 

C
p 


(9)
(6)
(10)
(8)
(11)
1
1
MATLAB Básico
Traduzindo as equações diferenciais para o Matlab:
Matlab
Real
dy(1)
dh/dt
y(1)
h
dy(2)
dT/dt
y(2)
T
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10;
% m3/h
Cp = 0.75; % kJ/(kg . K)
Ro = 1000; % kg/m3
U = 150;
% kJ/(m2 . s . K)
Te = 530;
%K
Th = 540;
%K
MATLAB Básico
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; % h
% Simulação do modelo
[t,y]=ode45('dydt',t,[(5/A) Th],[],[U A Ro Cp Fe R Te Th]);
MATLAB Básico
% Visualização da simulação
figure(1);
plot(t,y(:,1));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Altura (m)');
figure(2);
plot(t,y(:,2));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Temperatura (K)');
A única modificação em relação ao exemplo anterior é que
estamos passando duas condições iniciais (pois existem duas
variáveis dependentes):
MATLAB Básico
[t,y]=ode45('dydt',t,[(5/A) Th],[],[U A Ro Cp Fe R Te Th]);
MATLAB Básico
A função .m tem o código apresentado a seguir:
function dy = dydt(t,y,flag,par);
U = par(1);
A = par(2);
Ro = par(3);
Cp = par(4);
Fe = par(5);
R = par(6);
Te = par(7);
Th = par(8);
dy(1) = (Fe-(y(1)/R))/A;
dy(2) = (1/y(1))* ( ((Fe*Te/A)+(U*Th/(Ro*Cp)))...
- ( y(2)*((Fe/A)+(U/(Ro*Cp)))) );
dy = dy(:);
O vetor dy é criado como vetor linha (dy(1)) e
(dy(2)). Porém temos que retornar como vetor coluna.
Use o comando:
MATLAB Básico
matriz coluna = matriz linha (:)
Dica!
Quando for fazer os gráficos no programa principal
lembre-se que a primeira coluna de “dy” refere-se a “h” e a
segunda a “T”. Então para graficar h vs. tempo faça:
MATLAB Básico
figure(1);
plot(t,y(:,1));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Altura (m)');
Dica!
MATLAB Básico
Exemplo
4
Na compra de uma calculadora gráfica, a loja ofereceu duas propostas de
financiamento – proposta A e B. A proposta A é composta por 7 parcelas mensais
iguais de 114 reais cada. Já a proposta B prevê 10 parcelas mensais iguais de 98 reais
cada. Qual é a melhor opção de compra considerando a taxa de juros oferecida em
investimentos denominados “renda fixa”?
MATLAB Básico
A princípio poderia resolver o problema simplesmente multiplicado 114 x 7 e 10 x
98, achando o valor final pago. Os valores encontrados seriam 798 e 980. Logo, a
Proposta A parece mais favorável para o comprador.
É importante lembrar, porém, que essa forma de resolução não considera que o
dinheiro desvaloriza-se ao longo dos meses. Ou seja, o poder de compra de 100 reais
hoje, é superior ao poder de compra de 100 reais daqui a 10 meses. Outra forma de
pensar é considerar o “custo de oportunidade” – a taxa de retorno livre de que
conseguiria para o meu dinheiro caso, ao invés de pagar agora, investisse. De uma
forma ou de outra o que precisamos é do VALOR PRESENTE (VP) de cada série de
pagamentos, sendo os pagamentos descontados a dada taxa de juros. Para trazer
VALOR FUTURO (VF) para valor presente usa-se a fórmula:
VP = VF / ( 1 + i )n
Onde “i” é a taxa de juros mensal e “n” o número de meses entre o VF e o VP.
MATLAB Básico
clc
close all
clear all
ivetor=0:0.01:0.50;
MATLAB Básico
VPvetor114=[];
VPvetor98=[];
prompt{1}='Número de meses do pagamento da serie A:';
prompt{2}='Número de meses do pagamento da serie B:';
prompt{3}='Valor de cada parcela da serie A:';
prompt{4}='Valor de cada parcela da serie B:';
resposta=inputdlg(prompt,'Calculo da taxa de equilibrio');
nummeses114=str2num(char(resposta(1)));
nummeses98=str2num(char(resposta(2)));
v114=str2num(char(resposta(3)));
v98=str2num(char(resposta(4)));
MATLAB Básico
for J = 1:length(ivetor),
i=ivetor(J);
VP=[];
for K = 1:nummeses114,
VP(K)=v114/(1+i)^K;
end
VPfinal=sum(VP);
VPvetor114=[VPvetor114, VPfinal];
VP=[];
for K = 1:nummeses98,
VP(K)=v98/(1+i)^K;
end
VPfinal=sum(VP);
VPvetor98=[VPvetor98, VPfinal];
end
MATLAB Básico
plot(ivetor*100,VPvetor114,'-b')
hold on
plot(ivetor*100,VPvetor98,'-r')
title('Valor presente das parcelas a serem pagas')
legend( [ num2str(nummeses114), ' parc de ', num2str(v114) ,' reais cada'] , ...
[ num2str(nummeses98), ' parc de ', num2str(v98) ,' reais cada' ] )
xlabel('Taxa de juros mensal')
ylabel('Valor presente em Reais')
if (VPvetor114(1)<VPvetor98(1)),
posicoes=VPvetor114<VPvetor98;
achaZero=find(posicoes==0);
achaPrimeiroZero=min(achaZero);
plot(100*ivetor(achaPrimeiroZero),VPvetor114(achaPrimeiroZero),'ok')
else
posicoes=VPvetor98<VPvetor114;
achaZero=find(posicoes==0);
achaPrimeiroZero=min(achaZero);
plot(100*ivetor(achaPrimeiroZero),VPvetor114(achaPrimeiroZero),'ok')
end
text(100*ivetor(achaPrimeiroZero),50+VPvetor114(achaPrimeiroZero), ...
['Juros de equilibrio (a.m.) = ',num2str(100*ivetor(achaPrimeiroZero)),' %'] )
MATLAB Básico
Exemplo
5
Determinado processo possui função custo definida pela equação:
Y=((x-3)2)-6
MATLAB Básico
É necessário encontrar x que minimize o valor de Y. Não é difícil de visualizar
que a solução do problema é fazer x=3, de modo a levar Y um mínimo (-6). Mesmo
sabendo previamente a solução, vamos resolver através do MATLAB. Utilizamos
então a função “fminsearch”.
%Calculo do valor de x que minimiza a funcao custo
xmin = fminsearch('((x-3).^2)-6', 4)
%Gráfico da funcao custo
x=0:0.01:5;
y=((x-3).^2)-6;
plot(x,y);
hold on
MATLAB Básico
%Marca o ponto de minimo:
ymin=((xmin-3).^2)-6;
plot(xmin, ymin,'or')
MATLAB Básico
Ou...
MATLAB Básico
%Calculo do valor de x que minimiza a funcao custo
%Gráfico da funcao custo
x=0:0.01:5;
y=((x-3).^2)-6;
plot(x,y);
hold on
drawnow
function [y] = custo(x)
xmin = fminsearch('custo', 4)
plot(x, y,'ob')
hold on
pause(0.1)
%Marca o ponto de minimo:
ymin=((xmin-3).^2)-6;
plot(xmin, ymin,'or')
y=((x-3).^2)-6;
MATLAB Básico
Exemplo
6
Quando ajustamos uma curva a um conjunto de pontos experimentais, estamos minimizando a
distância entre a curva e os dados. Definindo essa distância como “erro”, estamos manipulando os
parâmetros que definem a curva de modo a minimizar o erro. Nesse caso “erro” é a minha “função
objetivo” a ser minimizada. É através dessa ótica que se torna possível usar “fminsearch” para
encontrar o valor ótimo dos parâmetros de ajuste de uma curva aos dados experimentais.
global yexp xexp
MATLAB Básico
%pontos experimentais
yexp=[1.1 2.12 2.85 4.4 5.0 6.5];
xexp=[1 2 3 4 5 6];
Parametros = fminsearch('custo',[1,2]);
function [saida] = custo(x)
global yexp xexp
a=x(1);
b=x(2);
yteo=a.*xexp + b; %calcula o valor teorico
%para cada pto experimental
MATLAB Básico
yerro=abs(yexp-yteo); %calculo do erro
saida=sum(yerro);
plot(xexp,yexp,'r*',xexp,yteo,'b-')
drawnow
pause(0.3)
MATLAB Básico
Exemplo
7
MATLAB Básico
MATLAB Básico
As equações diferenciais que descrevem
o processo são:
O modelo matemático do nosso reator CSTR tende ao estado estacionário.
Ou seja, seus parâmetros tendem a ficar constantes no tempo infinito.
Seria interessante introduzir perturbações em algumas variáveis e observar
como o reator se comporta.
Uma perturbação degrau em uma entrada u do sistema é tal
que:
MATLAB Básico
u = u0 ,
u = u0 + du,
t < tdegrau
t > tdegrau
Ou seja: antes do degrau a entrada u vale u0. Após o tempo
determinado para que o degrau ocorra (tdegrau) temos que u passa a
valer u0 + du.
Programa principal:
MATLAB Básico
% Definição das constantes do modelo
U =50; % BTU/(h.ft2.R)
A = 120; % ft2
DH = -30000; % BTU/lbm
Ro = 50; % lb/ft3
Cp = 0.75; % BTU/(lbm.R)
E = 30000; % BTU/lbm R = 1.99; % BTU/(lbm.R)
k0 = 7.08e10; % 1/h
V =48; % ft3
Te = 580; %R
Th = 550; %R
Fe = 18; % ft3/h
Cre = 0.48; % lbm/ft3
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; %h
% Perturbação na vazão de entrada
td = 5.0; %Tempo onde ocorre o degrau
fd = 2*Fe; %Valor assumido após o degrau
continua...
MATLAB Básico
Programa principal (continuação):
% Condições iniciais
Cr0 = 0.16; % lbm/ft3
T0 = 603; %R
% Simulação do modelo
[t,y] = ode45('dcstrdeg',t,[Cr0 T0],[],[U A DH Ro Cp E R k0 V Te Th …
Fe Cre],[td fd]);
% Visualização da simulação
figure(1);
plot(t,y(:,1)); title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)'); ylabel('Concentração de Reagente (lbm/ft3)');
figure(2);
plot(t,y(:,2)); title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)'); ylabel('Temperatura (R)');
Função “dcstrdeg”:
function dy = dcstrdeg(t,y,flag,par,deg);
U = par(1); A = par(2);
DH = par(3); Ro = par(4);
Cp = par(5); E = par(6);
R = par(7); k0 = par(8);
V = par(9); Te = par(10);
Th = par(11);
MATLAB Básico
continua...
Função “dcstrdeg” (continuação):
%Verifica a ocorrência de degrau:
if t >= deg(1)
Fe = deg(2);
else
Fe = par(12);
end;
Cre = par(13);
MATLAB Básico
dy(1) = (Fe/V)*(Cre-y(1)) - k0*exp(-E/(R*y(2)))*y(1);
dy(2) = (Fe/V)*(Te-y(2)) + ((DH*k0*exp(-E/(R*y(2)))*y(1))/(Ro*Cp)) - ...
(U*A*(y(2)-Th)/(V*Ro*Cp));
dy = dy(:);
MATLAB Básico
Exemplo
8
MATLAB Básico
Uma das grandes vantagens no uso de ferramentas computacionais é
reduzir o nosso esforço repetitivo, tarefa para a qual o computador é
muito eficiente. Supomos que temos um processo no qual
gostaríamos de testar uma série de condições iniciais. Para cada nova
condição inicial teríamos de refazer todas as contas. Um esforço
enorme! As linguagens de programação, e o Matlab em particular,
resolvem esse problema facilmente usando o já apresentado
comando “for”.
Usaremos o mesmo reator CSTR do exemplo anterior.
Programa principal:
MATLAB Básico
% Definição das constantes do modelo
U =50; % BTU/(h.ft2.R)
A = 120; % ft2
DH = -30000; % BTU/lbm
Ro = 50; % lb/ft3
Cp = 0.75; % BTU/(lbm.R)
E = 30000; % BTU/lbm R = 1.99; % BTU/(lbm.R)
k0 = 7.08e10; % 1/h V =48; % ft3
Te = 580; %R
Th = 550; %R
Fe = 18; % ft3/h
Cre = 0.48; % lbm/ft3
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; %h
% Condições iniciais
Cr0 = [0.16 0.32 0.48 0.64]; % lbm/ft3
T0 = 603; %R
continua...
Programa principal (continuação):
MATLAB Básico
% Simulação e visualização do modelo em batelada
cor = 'brmk';
leg = ['Cr0=0.16'; 'Cr0=0.32'; 'Cr0=0.48'; 'Cr0=0.64'];
for aux = 1 : length(Cr0)
[t,y] = ode45('dcstr',t,[Cr0(aux) T0],[], [U A DH Ro Cp E R k0 V…
Te Th Fe Cre]);
% Visualização da simulação
figure(1); hold on;
plot(t,y(:,1),cor(aux)); title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)');
ylabel('Concentração de Reagente (lbm/ft3)');
figure(2); hold on;
plot(t,y(:,2),cor(aux)); title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)');
ylabel('Temperatura (R)');
end;
continua...
Programa principal (continuação):
figure(1); legend(leg);
figure(2); legend(leg);
hold off;
MATLAB Básico
A seqüência de cores usadas é dada pelo vetor
“cor”, “letra a letra” através da flag. Consulte o
comando “plot” para detalhes.
Dica!
MATLAB Básico
Carlos André Vaz Junior
[email protected]
http://www.eq.ufrj.br/links/h2cin/carlosandre