Carlos Andr

1
Carlos André Vaz Junior cavazjunior_at_gmail.com http
//www.eq.ufrj.br/links/h2cin/carlosandre
2
O mundo MATLAB
3
Ajuda
?
4
Livros
5
Ambiente de Trabalho
6
Ambiente de Trabalho
7
Command Window
8
Command Window
Agora a 2, faço tudo de novo?!
9
Arquivo de Programação m-file
10
Current Directory
11
Current Directory
12
Criando variáveis
Tipos Básicos
Matriz
Estrutura
Char Array
Case Sensitive!
CaSe SeNsItIvE!
13
Tipo Matriz
Criando uma matriz
14
Tipo Char Array
Criando um char array
15
Tipo Estrutura
Banco de Dados da Turma Alunos Carla,
João, Bruno, Luis, Marcela Professor Marcelo
Horário 13h Sala 221
Estrutura
turma.alunos.nomesstrvcat( 'carla',joao','bruno'
, ... 'luis', 'marcela ) turma.professor.nome(
Marcelo) turma.horario1300 turma.sala221
16
Comando who e whos
Comando who e whos
17
Dicas!
Use para evitar que o resultado apareça na
tela.
Use A00.510 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.
Use size(A) para identificar as dimensões da
matriz. A maior dimensão é dada pelo comando
length(A)
Dicas!
18
Operações Matemáticas Simples
i) Soma e subtração soma (ou subtrai) elemento
por elemento da matriz. AB A-B   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
19
Operações Matemáticas Simples
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) 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)  
20
Operações Matemáticas Simples
ix) Determinante det(matriz)   x)
Inversa inv(matriz)   xi) Dimensões da
matriz size(matriz) lenght(matriz) numel(matriz)
 
Veja também flipud e fliplr
21
Referenciar um Elemento de uma Matriz
1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16
Elemento Matriz(2,3) ou Matriz(10)
22
Arquivo Function
23
Escopo das Variáveis
Função Alfa
Função Beta
A1 B2 global C C100
E15 F55 C23
24
Exemplo Rápido
Achando a posição do menor valor de uma matriz
x1 2 3 4 5 6 2 1 3 3 2 1 Forma
linear xminmin(x)
xminmin(xmin) i,jfind(xxmin)
Forma condensada i,jfind(x(min(min(x
))))
25
Exemplo Rápido
Achando o menor valor de uma função
X fzero('sin',2)
função
estimativa inicial
Veja também fsolve e fmin
26
Estruturas Lógicas
if
27
Estruturas Lógicas
AND
OR
28
Estruturas Lógicas
OR
a
b
resultado
1
1
1
0
1
1
1
1
0
0
0
0
29
Estruturas Lógicas
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
30
Estruturas Lógicas
  While   while I lt 10,
disp(oi) II1
end
Manipule o ponteiro I na rotina executada pelo
while
31
Estruturas Lógicas
     For   for J 1100,
A(1,J) 1/(IJ-1) end
Incremento automático do ponteiro J a cada loop
32
Proteção Contra Erros
     Try   try I 15
J teste A 1/(IJ-1) catch
disp(Erro na divisão) end
33
Encerrando uma Rotina
Break   i0 while i lt 100,
ii1 disp(i)
if igt10, break
end end
34
Gráficos
gtgt figure(1)
gtgt figure(2)
gtgt t00.0110 gtgt ysin(t) gtgt plot(t,y)
gtgt zcos(t) gtgt plot(t,z)
35
Gráficos
gtgt figure(3)
gtgt plot(t,y)
gtgt subplot(1,2,2) gtgt plot(t,z)
gtgt subplot(1,2,1)
36
Gráficos
gtgt t00.2510 gtgt ysin(t) gtgt plot(t,y,'r') gtgt
xlabel('tempo') gtgt ylabel('seno') gtgt title('Seno
vs. Tempo') gtgt Axis(0 10 -2 2)
37
Gráficos
gtgt t00.0110 gtgt ysin(t) gtgt zcos(t) gtgt
plot(t,y,'g-',t,z,'r-') gtgt legend('seno','cosseno'
)
Ou...
gtgt t00.0110 gtgt ysin(t) gtgt zcos(t) gtgt
plot(t,y,'g-) gtgt hold on gtgt plot(t,z,'r-') gtgt
legend('seno','cosseno')
38
Gráficos - Tortas
gtgt x 1 3 0.5 2.5 2 gtgt explode 0 1 0 0
0 gtgt pie(x,explode) gtgt colormap jet gtgt
legend('EMB','IND','ACO','DIV','POT')
39
Gráficos - Tortas
gtgt x 1 3 0.5 2.5 2 gtgt explode 0 1 0 0
0 gtgt pie3(x,explode) gtgt colormap jet gtgt
legend('EMB','IND','ACO','DIV','POT')
40
Gráficos - Barras
gtgt x -2.90.22.9 gtgt bar(x,exp(-x.x)) gtgt
colormap hsv
41
Gráficos - Superfície
42
Gráficos - Superfície
xx00.011 yy00.011 X,Ymeshgrid(xx,yy)
Zexp(-0.5(X.2Y.2)) colormap
jet figure(1)surf(X,Y,Z) rotate3d on shading
interp
43
Gráficos - Superfície
44
Gráficos - Superfície
45
Gráficos - Superfície
xx00.011 yy00.011 X,Ymeshgrid(xx,yy)
ZX.X.Y colormap jet figure(1)surf(X,Y,Z)
rotate3d on shading interp
46
Gráficos - Superfície
47
Gráficos - Superfície
48
Gráficos - Superfície
Malha triangular da base malha da base
xx00.011 yy00.011 X,Ymeshgrid(xx,yy)
Z1-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")
izfind(Zlt0)Z(iz)nan
49
Gráficos - Superfície
Composição (3 componentes)
Malha triangular da base malha da base
xx00.011 yy00.011 X,Ymeshgrid(xx,yy)
Z1-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")
izfind(Zlt0)Z(iz)nan
50
Gráficos - Superfície
Alguns Z são negativos! Não pode!
Malha triangular da base malha da base
xx00.011 yy00.011 X,Ymeshgrid(xx,yy)
Z1-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")
izfind(Zlt0)Z(iz)nan
51
Gráficos - Superfície
vv1(X.log(X))(Y.log(Y))(Z.log(Z))
gráfico da superfície colormap
jet figure(1)surf(X,Y,vv1) rotate3d on shading
interp xlabel('X1')ylabel('X2')zlabel('DeltaGi/
RT')
52
Gráficos - Superfície
53
Gráficos
Use x,yginput(2) para capturar dois pontos
no gráfico
Use close all para fechar todas as figuras
Use clf para apagar a figura atual
Dica!
54
Exemplos
55
Exemplo 1
56
Modelagem simples de um tanque de nível
57
Modelagem simples de um tanque de nível
58
Modelagem simples de um tanque de nível
59
Modelagem simples de um tanque de nível
60
Modelagem simples de um tanque de nível
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
RFe(1 - exp(-t/(RA))) m
Visualização da simulação plot(t,h) title('Simula
ção do tanque de nível') xlabel('Tempo
(h)') ylabel('Altura (m)')
61
Modelagem simples de um tanque de nível
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!
62
Exemplo 2
63
Modelagem de um tanque de nível via ED
Muitas vezes é muito trabalhoso, ou mesmo
impossível, encontrar a solução analítica para o
conjunto de equações diferenciais. Nesse caso
temos que simular usando solução numérica das
equações diferenciais. Vamos assumir que o
modelo do exemplo 1 não tivesse solução
analítica, e então usar o Matlab para estudar o
comportamento da altura do nível com o tempo. A
equação diferencial será
64
Modelagem de um tanque de nível via ED
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
65
Modelagem de um tanque de nível via ED
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 EDOs por método de
Runge-Kutta de ordem média (consulte o help sobre
ODE45 para maiores detalhes). t,h
ode45('dhdt',t, 0,,R A Fe)
66
Modelagem de um tanque de nível via ED
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. -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 EDOs. 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.
67
Modelagem de um tanque de nível via ED
-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.   -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).
68
Modelagem de um tanque de nível via ED

• 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 EDOs (todas as
  ODExx), a declaração deste arquivo deve seguir a
  sintaxe
• 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.

69
Exemplo 3
70
Modelagem de um tanque de aquecimento
71
Modelagem de um tanque de aquecimento
72
Modelagem de um tanque de aquecimento
Traduzindo as equações diferenciais para o Matlab
73
Modelagem de um tanque de aquecimento
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   Tempo de
simulação t 0.0 0.01 10.0 h  
Simulação do modelo t,yode45('dydt',t,(5/A)
Th,,U A Ro Cp Fe R Te Th)  
74
Modelagem de um tanque de aquecimento
  Visualização da simulação figure(1) plot(t,y(
,1)) title('Tanque de aquecimento') xlabel('Tem
po (h)') ylabel('Altura (m)') figure(2) plot(t
,y(,2)) title('Tanque de aquecimento') xlabel(
'Tempo (h)') ylabel('Temperatura (K)')
75
Modelagem de um tanque de aquecimento
A única modificação em relação ao exemplo
anterior é que estamos passando duas condições
iniciais (pois existem duas variáveis
dependentes) t,yode45('dydt',t,(5/A)
Th,,U A Ro Cp Fe R Te Th)
76
Modelagem de um tanque de aquecimento
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)) (
((FeTe/A)(UTh/(RoCp)))... - (
y(2)((Fe/A)(U/(RoCp)))) ) dy dy()
77
Modelagem de um tanque de aquecimento
O vetor dy é criado como vetor linha (dy(1)) e
(dy(2)). Porém temos que retornar como vetor
coluna. Use o comando matriz coluna
matriz linha ()
Dica!
78
Modelagem de um tanque de aquecimento
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
figure(1) plot(t,y(,1)) title('Tanque de
aquecimento') xlabel('Tempo (h)')
ylabel('Altura (m)')
Dica!
79
Exemplo 4
80
Aplicando perturbações no CSTR
81
Aplicando perturbações no CSTR
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.
82
Aplicando perturbações no CSTR
Uma perturbação degrau em uma entrada u do
sistema é tal que   u u0 , t lt tdegrau u
u0 du, t gt 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.
83
Aplicando perturbações no CSTR
Programa principal
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 2Fe
Valor assumido após o degrau
continua...
84
Aplicando perturbações no CSTR
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)')
85
Aplicando perturbações no CSTR
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)  
continua...
86
Aplicando perturbações no CSTR
Função dcstrdeg (continuação)
Verifica a ocorrência de degrau if t gt deg(1)
Fe deg(2) else Fe
par(12) end  Cre par(13)  dy(1)
(Fe/V)(Cre-y(1)) - k0exp(-E/(Ry(2)))y(1) dy(2
) (Fe/V)(Te-y(2)) ((DHk0exp(-E/(Ry(2)))y(
1))/(RoCp)) - ... (UA(y(2)-Th)/(VR
oCp)) dy dy()
87
Exemplo 5
88
Simulação em batelada
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.
89
Simulação em batelada
Programa principal
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...
90
Simulação em batelada
Programa principal (continuação)
Simulação e visualização do modelo em
batelada cor 'brmk' leg 'Cr00.16'
'Cr00.32' 'Cr00.48' 'Cr00.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...
91
Simulação em batelada
Programa principal (continuação)
figure(1) legend(leg) figure(2)
legend(leg) hold off
A seqüência de cores usadas é dada pelo vetor
cor, letra a letra através da flag. Consulte
o comando plot para detalhes.
Dica!
92
Carlos André Vaz Junior cavazjunior_at_gmail.com http
//www.eq.ufrj.br/links/h2cin/carlosandre