Embed Size (px)
Citation preview
MA
TLA
B B
ásic
o
Mais de 1 milhão de resultados
MA
TLA
B B
ásic
o
?http://newsreader.mathworks.com
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
Agora a = 2, faço tudo de novo?!
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
Char ArrayMatriz
Tipos Básicos
Case Sensitive!
Estrutura
CaSe SeNsItIvE!
MA
TLA
B B
ásic
o
Criando uma matriz:
MA
TLA
B B
ásic
o
Criando um “char array”:
MA
TLA
B B
ásic
o
Estrutura:turma.alunos.nomes=strvcat( 'carla',’joao','bruno', ... 'luis', 'marcela‘ );turma.professor.nome=(‘Marcelo‘)turma.horario=1300turma.sala=221
Banco de Dados da “Turma”: Alunos: Carla, João, Bruno, Luis, Marcela Professor: Marcelo Horário: 13h Sala: 221
MA
TLA
B B
ásic
o
Comando “who” e “whos”
MA
TLA
B B
ásic
o
Use A=0:0.5:10 para gerar matrizes com dados em seqüência.
Use “;” para evitar que o resultado apareça na tela.
Dicas!
Use “clear A” para apagar a variável A.
Use “size(A) ” para identificar as dimensões da matriz. A maior dimensão é dada pelo comando “length(A) ”
Use “clear all” para apagar todas as variáveis armazenadas.
MA
TLA
B B
ásic
o
i) Soma e subtração: soma (ou subtrai) elemento por elemento da matriz. A+B 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 .* BA ./ B
MA
TLA
B B
ásic
o
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)
MA
TLA
B B
ásic
o
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
MA
TLA
B B
ásic
o
1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16
Elemento = Matriz(2,3) ou Matriz(10)
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
A=1B=2global CC=100
Função AlfaE=15F=55C=23
Função Beta
Programa Principal / Workspaceglobal CC=100D=22
MA
TLA
B B
ásic
o
x=[1 2 3 4 5 6; 2 1 3 3 2 1];
%Forma linear: xmin=min(x); xmin=min(xmin); [i,j]=find(x==xmin);
Achando a posição do menor valor de uma matriz:
%Forma condensada: [i,j]=find(x==(min(min(x))));
MA
TLA
B B
ásic
o
X = fzero('sin',2)
Achando o menor valor de uma função:
função estimativa inicial
Veja também: fsolve e fmin
MA
TLA
B B
ásic
o
if:
MA
TLA
B B
ásic
o
AND OR
MA
TLA
B B
ásic
o
Falso Verdadeiro
AND
a b resultado1 1 1
0 1 0
1 0 0
0 0 0
OR
a b resultado1 1 1
0 1 1
1 0 1
0 0 0
MA
TLA
B B
ásic
o
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
MA
TLA
B B
ásic
o
While: while I < 10, disp(‘oi’); I=I+1; end
Manipule o ponteiro I na rotina executadapelo “while”
MA
TLA
B B
ásic
o
For: for J = 1:100, A(1,J) = 1/(I+J-1); end
Incremento automático do ponteiro J a cada loop
MA
TLA
B B
ásic
o
>> figure(1)
>> t=0:0.01:10;>> y=sin(t);>> plot(t,y)
>> figure(2)
>> z=cos(t);>> plot(t,z)
MA
TLA
B B
ásic
o
Use “close all” para fechar todas as figuras
Dica!
Use “clf” para apagar a figura atual
Use “[x,y]=ginput(2)” para capturar dois pontos no gráfico
MA
TLA
B B
ásic
o
>> figure(3)
>> subplot(1,2,1)
>> plot(t,y)
>> subplot(1,2,2)>> plot(t,z)
MA
TLA
B B
ásic
o
>> 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])
MA
TLA
B B
ásic
o
>> t=0:0.01:10;>> y=sin(t);>> z=cos(t);>> plot(t,y,'g-',t,z,'r-')>> legend('seno','cosseno')
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')
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
xx=0:0.01:1;yy=0:0.01:1;[X,Y]=meshgrid(xx,yy);
Z=exp(-0.5*(X.^2+Y.^2));
colormap jetfigure(1);surf(X,Y,Z); rotate3d on; shading interp;
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
%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;
MA
TLA
B B
ásic
o
%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;
Composição(3 componentes)
MA
TLA
B B
ásic
o
%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ãonegativos! Não pode!
MA
TLA
B B
ásic
o
vv1=(X.*log(X))+(Y.*log(Y))+(Z.*log(Z));
%gráfico da superfíciecolormap jetfigure(1);surf(X,Y,vv1); rotate3d on; shading interp;xlabel('X1');ylabel('X2');zlabel('DeltaGi/RT');
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
Exemplos
MA
TLA
B B
ásic
o
Exemplo1
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
E x e m p l o 1 :M o d e l o s s i m p l e s - o t a n q u e d e n í v e l
1
h
F E
F
A
h
F E
F
A
h
F E
F
A
C o n s i d e r a n d o c o n s t a n t e s a v a z ã o d e a l i m e n t a ç ã o F E , a d e n s i d a d e e a t e m p e r a t u r a T , e q u e o s i s t e m a e s t á s u j e i t o àc o n d i ç ã o i n i c i a l :
00 hth ( 1 )
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - o t a n q u e d e n í v e lp o d e - s e e s c r e v e r o b a l a n ç o d e m a s s a d o s i s t e m a
1
FFdt
tdmE
dt
tdhAdt
tdm
FFAdt
tdhE
1
A i n d a ,
e , p o r t a n t o ,
( 2 )
( 3 )
( 4 )
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - o t a n q u e d e n í v e lF r e q ü e n t e m e n t e , c o n s i d e r a - s e a v a z ã o d e s a í d a d o t a n q u e
p r o p o r c i o n a l à a l t u r a d a c o l u n a d e l í q u i d o é i n v e r s a m e n t e p r o p o r c i o n a l a u m a r e s i s t ê n c i a a o e s c o a m e n t o ( R ) :
1
RhF
RhF
Adttdh
E1
L o g o ,
( 5 )
( 6 )
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - o t a n q u e d e n í v e lE s t e m o d e l o s i m p l e s d e u m t a n q u e d e n í v e l , s e m b a l a n ç o d e
e n e r g i a , p o s s u i u m a s o l u ç ã o a n a l í t i c a :
1
RA
t
E eRFth 1
P a r a s i m u l a r e s t e m o d e l o , b a s t a e s c o l h e r o s v a l o r e s d a s c o n s t a n t e s R , A e F E , d a s c o n d i ç õ e s i n i c i a i s h 0 e t 0 .
A s i m u l a ç ã o d a s o l u ç ã o a n a l í t i c a d o m o d e l o d o t a n q u e d e n í v e l é m o s t r a d a a s e g u i r .
( 7 )
MA
TLA
B B
ásic
o
% Definição das constantes do modeloR = 1; % h/m2A = 2; % m2Fe = 10; % m3/h% Tempo de simulaçãot = 0.0 : 0.01 : 10.0; % h% Simulação da altura de líquidoh = R*Fe*(1 - exp(-t/(R*A))); % m% Visualização da simulaçãoplot(t,h);title('Simulação do tanque de nível');xlabel('Tempo (h)');ylabel('Altura (m)');
MA
TLA
B B
ásic
o
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!
MA
TLA
B B
ásic
o
Exemplo2
MA
TLA
B B
ásic
o
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á:
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - o t a n q u e d e n í v e lF r e q ü e n t e m e n t e , c o n s i d e r a - s e a v a z ã o d e s a í d a d o t a n q u e
p r o p o r c i o n a l à a l t u r a d a c o l u n a d e l í q u i d o é i n v e r s a m e n t e p r o p o r c i o n a l a u m a r e s i s t ê n c i a a o e s c o a m e n t o ( R ) :
1
RhF
RhF
Adttdh
E1
L o g o ,
( 5 )
( 6 )
MA
TLA
B B
ásic
o
% Definição das constantes do modeloR = 1; % h/m2A = 2; % m2Fe = 10; % m3/h% Tempo de simulaçãot = 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çãoplot(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;
MA
TLA
B B
ásic
o
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).
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
MA
TLA
B B
ásic
o
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 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.
MA
TLA
B B
ásic
o
-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).
MA
TLA
B B
ásic
o
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:
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.
MA
TLA
B B
ásic
o
Exemplo3
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
E x e m p l o 3 T a n q u e d e a q u e c i m e n t o
1
h
F E , T E
F , T
A
T h
T h
h
F E , T E
F , T
A
T h
T h
C o n s i d e r a n d o c o n s t a n t e s a v a z ã o d e a l i m e n t a ç ã o F E , a t e m p e r a t u r a T h , o c o e fi c i e n t e g l o b a l d e t r a n s f e r ê n c i a d e c a l o r U e a s p r o p r i e d a d e s d o fl u i d o e C p e q u e o s i s t e m a e s t á s u j e i t o à s c o n d i ç õ e s i n i c i a i s :
00 hth 00 TtT
MA
TLA
B B
ásic
o
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - t a n q u e d e a q u e c i m e n t o
C o m o n o c a s o a n t e r i o r , o b a l a n ç o d e m a s s a p o d e s e r e s c r i t o c o m o :
1
RhF
Adttdh
E1
( 6 )
O b a l a n ç o d e e n e r g i a é e s c r i t o c o m o :
QFHHFdtVTdC EEp ( 8 )
M o d e l a g e m & D i n â m i c a d e P r o c e s s o s
M o d e l o s s i m p l e s - t a n q u e d e a q u e c i m e n t o
1
dtdhT
dtdThA
dtdVT
dtdTV
dtVTd
( 9 )
QFHHFRhF
AT
dtdThAC EEEp
( 1 0 )
p
E
p
hEE
CU
AFT
CUT
ATF
hdtdT
1
( 1 1 )
MA
TLA
B B
ásic
o
Matlab Real
dy(1) dh/dt
y(1) h
dy(2) dT/dt
y(2) T
Traduzindo as equações diferenciais para o Matlab:
MA
TLA
B B
ásic
o
% Definição das constantes do modeloR = 1; % h/m2A = 2; % m2Fe = 10; % m3/hCp = 0.75; % kJ/(kg . K)Ro = 1000; % kg/m3U = 150; % kJ/(m2 . s . K)Te = 530; % K Th = 540; % K % Tempo de simulaçãot = 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]);
MA
TLA
B B
ásic
o
% Visualização da simulaçãofigure(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)');
MA
TLA
B B
ásic
o
A única modificação em relação ao exemplo anterior é que estamos passando duas condições iniciais (pois
existem duas variáveis dependentes):
[t,y]=ode45('dydt',t,[(5/A) Th],[],[U A Ro Cp Fe R Te Th]);
MA
TLA
B B
ásic
o
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(:);
MA
TLA
B B
ásic
o
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!
MA
TLA
B B
ásic
o
figure(1);plot(t,y(:,1));title('Tanque de aquecimento');xlabel('Tempo (h)'); ylabel('Altura (m)');
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:
Dica!
MA
TLA
B B
ásic
o
Exemplo4
MA
TLA
B B
ásic
o
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”?
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.
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
clcclose allclear all ivetor=0:0.01:0.50; 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)));
MA
TLA
B B
ásic
o
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
MA
TLA
B B
ásic
o
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)),' %'] )
MA
TLA
B B
ásic
o
Exemplo5
MA
TLA
B B
ásic
o
Determinado processo possui função custo definida pela equação:
Y=((x-3)2)-6
É 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”.
MA
TLA
B B
ásic
o
%Calculo do valor de x que minimiza a funcao custoxmin = fminsearch('((x-3).^2)-6', 4) %Gráfico da funcao custox=0:0.01:5;y=((x-3).^2)-6;plot(x,y);hold on %Marca o ponto de minimo:ymin=((xmin-3).^2)-6;plot(xmin, ymin,'or')
MA
TLA
B B
ásic
o
Ou...
MA
TLA
B B
ásic
o
%Calculo do valor de x que minimiza a funcao custo %Gráfico da funcao custox=0:0.01:5;y=((x-3).^2)-6;plot(x,y);hold ondrawnow xmin = fminsearch('custo', 4) %Marca o ponto de minimo:ymin=((xmin-3).^2)-6;plot(xmin, ymin,'or')
function [y] = custo(x) y=((x-3).^2)-6; plot(x, y,'ob')hold onpause(0.1)
MA
TLA
B B
ásic
o
Exemplo6
MA
TLA
B B
ásic
o
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 %pontos experimentaisyexp=[1.1 2.12 2.85 4.4 5.0 6.5];xexp=[1 2 3 4 5 6]; Parametros = fminsearch('custo',[1,2]);
MA
TLA
B B
ásic
o
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 yerro=abs(yexp-yteo); %calculo do errosaida=sum(yerro); plot(xexp,yexp,'r*',xexp,yteo,'b-')drawnowpause(0.3)
MA
TLA
B B
ásic
o
Exemplo7
MA
TLA
B B
ásic
o
MA
TLA
B B
ásic
o
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.
MA
TLA
B B
ásic
o
Uma perturbação degrau em uma entrada u do sistema é tal que:
u = u0 , t < tdegrauu = u0 + du, 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.
MA
TLA
B B
ásic
o
% Definição das constantes do modeloU =50; % BTU/(h.ft2.R)A = 120; % ft2DH = -30000; % BTU/lbmRo = 50; % lb/ft3Cp = 0.75; % BTU/(lbm.R)E = 30000; % BTU/lbm R = 1.99; % BTU/(lbm.R)k0 = 7.08e10; % 1/h V =48; % ft3Te = 580; %R Th = 550; %RFe = 18; % ft3/h Cre = 0.48; % lbm/ft3 % Tempo de simulaçãot = 0.0 : 0.01 : 10.0; %h % Perturbação na vazão de entradatd = 5.0; %Tempo onde ocorre o degraufd = 2*Fe; %Valor assumido após o degrau
Programa principal:
continua...
MA
TLA
B B
ásic
o
% Condições iniciaisCr0 = 0.16; % lbm/ft3T0 = 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çãofigure(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)');
Programa principal (continuação):
MA
TLA
B B
ásic
o
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);
Função “dcstrdeg”:
continua...
MA
TLA
B B
ásic
o
%Verifica a ocorrência de degrau:
if t >= deg(1)
Fe = deg(2);
else
Fe = par(12);
end;
Cre = par(13);
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(:);
Função “dcstrdeg” (continuação):
MA
TLA
B B
ásic
o
Exemplo8
MA
TLA
B B
ásic
o
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.
MA
TLA
B B
ásic
o
% Definição das constantes do modeloU =50; % BTU/(h.ft2.R)A = 120; % ft2DH = -30000; % BTU/lbmRo = 50; % lb/ft3Cp = 0.75; % BTU/(lbm.R)E = 30000; % BTU/lbm R = 1.99; % BTU/(lbm.R)k0 = 7.08e10; % 1/h V =48; % ft3Te = 580; %R Th = 550; %RFe = 18; % ft3/h Cre = 0.48; % lbm/ft3 % Tempo de simulaçãot = 0.0 : 0.01 : 10.0; %h % Condições iniciaisCr0 = [0.16 0.32 0.48 0.64]; % lbm/ft3T0 = 603; %R
Programa principal:
continua...
MA
TLA
B B
ásic
o
% Simulação e visualização do modelo em bateladacor = '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;
Programa principal (continuação):
continua...
MA
TLA
B B
ásic
o
figure(1); legend(leg);figure(2); legend(leg);hold off;
Programa principal (continuação):
A seqüência de cores usadas é dada pelo vetor “cor”, “letra a letra” através da flag. Consulte o comando “plot” para detalhes.
Dica!
