C03 SEANL v04

Prévia do material em texto

ENG-D01: Métodos 
Computacionais na Engenharia
Profa Karen PontesProf . Karen Pontes
Dezembro 2008
Métodos Computacionais na EngenhariaMétodos Computacionais na Engenharia
Raízes de Equações não linearesRaízes de Equações não-lineares
x0,f(x0)0 1l AP B
x
x1,f(x1)
x
0 ( 0)
raiz
0 1
1 1lnP BT
 
3 2 2( ) 0Z Z A B B Z AB  x0x1x2y=f(x
)
( ) 0Z Z A B B Z AB     
Raízes de Equações não linearesRaízes de Equações não-lineares
E ilíb i d f 0 1A Equilíbrio de fases: 0 11 1ln
AP B
T
 
RT Equação de estado:
( )
RT aP
V b V V b
  
3 2 2( ) 0Z Z A B B Z AB     
 Constante de reação:
( ) 0Z Z A B B Z AB
/E RTr Ae Constante de reação:
F t d f i ã d D i t b l t
r Ae
 Fator de fricção de Darcy para o regime turbulento:
1 / 2.51( ) 0 86858ln DF f    0.5 0.5
Re
( ) 0.86858ln
3.7D D D
F f
f N f
    
Raízes de Equações ReaisRaízes de Equações Reais
E ã l éb i ã li ( ) 0f Equação algébrica não linear:
 Objetivo: calcular x (raiz ou zero)
( ) 0f x 
 Objetivo: calcular x (raiz ou zero)
 Exemplo: exp( / ) 0x x c  p( )
Solução analítica: 1
1.1
analítica
é i
    
 
1 ln 1 1x c u c u u         0.8
0.9
c
 
v
a
l
u
e
numérica
 ln 1 1u c 
0 5
0.6
0.7
0.4 0.45 0.5 0.55 0.6
0.5
root x value
Raízes de Equações não linearesRaízes de Equações não-lineares
 Algumas vezes não é possível obter a solução analítica.
 Se possível, requer tempo e habilidade!
Métodos Pode-se desejar maior precisão. Métodos numéricos
 Raízes coincidentes
100
1.5
3 2( 1) ( 2) ( 3) 0x x x    exp( /10)sen(10 ) 0x x 
50
100
0.5
1
0
f
(
x
)
-0.5
0
f
(
x
)
-2 0 2 4
-50
x
5 10 15 20
-1.5
-1
x
Métodos GráficosMétodos Gráficos
N f l õ di t t Nem sempre fornecem soluções diretamente
 Podem fornecer estimativa inicial para o método iterativo 
Matlab: plot+ginput(n )Matlab: plot+ginput(nraízes)
( ) 0x xf x e x e x     
1
1.5
xe
x
4
6
0.5
1
f
(
x
)
x
0 (0.56,0.57)P 
0
2
f
(
x
)
0 1 2 3
0
x
-1 0 1 2 3 4 5 6
-2
0
x
1 2 30,0565; 4,3468; 4,6694x x x   
Métodos numéricos interativos EANLMétodos numéricos interativos - EANL
Métodos de confinamentoMétodos de confinamento Métodos abertosMétodos abertos
a) b)
Algoritmo para métodos iterativosAlgoritmo para métodos iterativos
E ti tiEstimativa 
inicial x0
Avaliação
de f
Estimativa
Converge?
de xi
Sim Não
Converge?Solução
Chave doChave do 
problema !
Critérios de convergênciaCritérios de convergência
E l ti x x Erro relativo:
Mí i d f
1
1
i i
i
x x
x
 
 Mínimo de f: 2( )if x 
 Diferença: 1 3i ix x  
 Limite de interações (h): i h onde i é a interação
 Apenas um critério pode não garantir a convergência
Método da BisseçãoMétodo da Bisseção
S f( ) é tí i t l [ b] ( ) ( ) 0f f b Se f(x) é contínua no intervalo [a,b] e
 f(x) muda de sinal em [a,b], portanto há pelo menos 
uma raiz neste intervalo f(x)
( ) ( ) 0f a f b 
uma raiz neste intervalo, f(x)
Redução sucessiva do intervalo (fator 1/2) Redução sucessiva do intervalo (fator 1/2)
Método da BisseçãoMétodo da Bisseção
1) ( ) ( ) 0f a f b  1 21 2
a bp  
p2
raiz
a3=
) ( ) ( )f f
1 12) ( ) ( ) 0
( ) ( ) 0
f p f b
f f
 

1 2
p
raiz  b
b=b1
a=a1
p1=b3
a2=
b2 =
3
1 1 ( ) ( ) 0f a f p  raiz  2 1
2 1
b p
a a


2 2
23) 2
a bp 
a1 b1
a b
2 2
2 2
( ) ( ) 0
( ) ( ) 0
f a f p
f p f b
 
 
2
raiz  3 2a pa2 b2
a3 b3
2 2( ) ( )f p f 3 2
3 2
p
b b
Método da BisseçãoMétodo da Bisseção
E l 1)(0 f)5/()(f Exemplo
3297,1)(;2
1)(;0


bfb
afa)5/exp()( xxxf 
]1;0[ 18127.0)1( 1  fp
]1;75,0[ 11071.0)75.0( 75.0
]1;5,0[ 40484.0)5.0( 5.0


fp
fp
]875,0;8125,0[ 037516.0)8125.0( 8125.0
]975,0;75,0[ 035543.0)5.0( 875.0


fp
fp
0,038e 107,9)84375.0( 84375.0
][)(
4  fp
fp
Método da BisseçãoMétodo da Bisseção
S t d d Sempre converge para uma resposta, desde que uma 
raiz esteja contida no intervalo inicial
 Problemas numéricos de divisão por zero nao ocorrem
 Tamanho do passo pode ser ajustadoTamanho do passo pode ser ajustado 
 Nem sempre o intervalo é conhecido a priori
 Método muito lento em comparação com outros 
métodos
Algoritmo para método da bisseção
clear all
F = inline('8-4.5*(x-sin(x))');
a = 2; b = 3; imax = 20; tol = 0.001;
Fa=F(a); Fb=F(b);Fa F(a); Fb F(b);
if Fa*Fb > 0
disp('Erro: A funcao tem o mesmo sinal nos pontos a e b.')
else
disp('iteracao a b (xNS) Solucao f(xNS) Tolerancia’)
for i = 1:imaxfor i = 1:imax
xNS = (a + b)/2;
toli=(b-a)/2;
FxNS=F(xNS);
fprintf('%3i %11.6f %11.6f %11.6f %11.6f %11.6f\n',i, a, b, xNS, FxNS, toli)
if FxNS == 0
fprintf(' Solucao exata x =%11.6f encontrada',xNS)
break
end
if toli < tolif toli tol
break
end
if i == imax
fprintf('Solucao nao foi obtida em %i iteracoes',imax)
breakbreak
end
if F(a)*FxNS < 0
b = xNS;
else
a = xNS;
end
end
end
Método da Posição Falsa ou Regula Falsi
 Dado [a b] f(x) é contínua e a equação possui solução
Método da Posição Falsa ou Regula Falsi
 Dado [a, b], f(x) é contínua e a equação possui solução
 Se  f(x) muda de sinal há solução( ) ( ) 0f a f b  Se  f(x) muda de sinal  há solução( ) ( ) 0f a f b 
Método da Posição Falsa ou Regula FalsiMétodo da Posição Falsa ou Regula Falsi
Cál l d i i d Cálculo da nova estimativa xNS: equação da reta que 
une os extremos a e b do intervalo
( ) ( ) ( ) ( )f b f a b f b ( ) ( ) ( ) ( )f fy x b f b
b a
   Reta passa por (a,f(a)) e (b,f(b))
(xNS , 0)
( ) ( )
NS
af b bf ax  0 1 1 02 ( ) ( )x f x x f xx 
0 1 2, e NSa x b x x x  
( ) ( )NS f b f a 2 1 0( ) ( )
x
f x f x
Método da Posição Falsa ou Regula FalsiMétodo da Posição Falsa ou Regula Falsi
0 1 1 0( ) ( )x f x x f x  ( )f x x x0 1 1 0
2
1 0
( ) ( )
( ) ( )
x f x x f xx
f x f x
 
 1 1 0
2 1
1 0
( )
( ) ( )
f x x x
x x
f x f x
  
 Generalizando:
1
1
1
( )( )
( ) ( )
n n n
n n
n n
f x x xx x
f x f x



  
Método da Posição FalsaMétodo da Posição Falsa
E l 1)(0 f)5/()(f Exemplo
3297,1)(;2
1)(;0
11
00


xfx
xfx)5/exp()( xxxf 
24
2
2
106,1e 101,2)84476,0( 84476,0
]85949,0;0[ 016253.0)85949,0( 85949,0
 

fx
fx
Método da Posição Falsa ou Regula FalsiMétodo da Posição Falsa ou Regula Falsi
C bi ã d ét d Combinação dos métodos:
 Bisseção: se , então a raiz está n intervalo
[a b]
( ) ( ) 0f a f b 
[a,b]
 Newton-Raphson: fórmula recursiva sem necessidade Newton Raphson: fórmula recursiva sem necessidade 
do cálculo da derivada
 Sempre converge
 Função côncava ou convexa: um dos extremos é o 
mesmo em todas as iterações. Solução poderia ser mais ç ç p
rápida
Método de Newton RaphsonMétodo de Newton-Raphson
Mét d i l Método mais popular
 Pode ser mais complicado  derivada, condição inicial
B d d i d f t l ã Baseado na derivada para fazer extrapolação
( ) 0 ( )f x f x
x0,f(x0)
0 0
0 1 0
0 1 0
( ) 0 ( )'( )
'( )
f x f xf x x x
x x f x
   
x1,f(x1)
raiz
1 1
1 2 1
1 2 1
( ) 0 ( )'( )
'( )
f x f xf x x x
x x f x
   
x0x1x2y=f(x) 1
( )
'( )
n
n n
nf xx x
f x
 
( )nf
Método de Newton RaphsonMétodo de Newton-Raphson
f( ) é d dif iá l d i d ã f(x) é duas vezes diferenciável, e suas derivadas são 
contínuas em [a,b]
 Expansão em série de Taylor de f(x) em torno da 
estimativa inicial x forneceestimativa inicial xn fornece
2
1
1 1
( )( ) ( ) ( ) '( ) ''( ) ...
2
n n
n n n n n n
x xf x f x x x f x f x 
      
 Se x é a solução f(x ) = 0 e x é um ponto próximo
2
 Se xn+1 é a solução, f(xn+1) = 0, e xn é um ponto próximo. 
Truncando a série no segundo termo:
( )f
10 ( ) ( ) '( )n n n nf x x x f x    1 ( )'( )
n
n n
n
f xx x
f x
 
Método de Newton RaphsonMétodo de Newton-Raphson
C ét d d N t di Casos em que o método de Newton diverge
Método de Newton RaphsonMétodo de Newton-Raphson
R ál l d d i d líti Requer o cálculo da derivada analítica 
D i d i d l l d i t Derivada aproximada calculada numericamente 
( ) ( )'( ) f x x f xf x   ( )f x
x
 
 Este método é muito bom, converge rapidamente, mas 
para problemas que possam convergir
Método de Newton RaphsonMétodo de Newton-Raphson
 Pode-se mostrar que o método converge quando:
 f(x), f ‘(x) e f “(x) forem contínuas
 f ‘(x) for diferente de zero na solução
 a estimativa inicial estiver próximo à solução exata
 Problemas de convergência ocorrem tipicamente Problemas de convergência ocorrem tipicamente 
quando o valor da derivada da função é próximo de zero 
na vizinhança da solução
Método de Newton RaphsonMétodo de Newton-Raphson
E l )5/()(f Exemplo
)5/exp(2,01)(';0 00 xxfx 
)5/exp()( xxxf 
01353,0 84458,0 1693,1)(' 013145,0)( 8333,0
8333,0 2,1)(' 1)( 0
6
2111
1000


exxfxfx
xxfxfx
0 84458,0 1689,1)(' 1014,2)( 8446,0 32
6
21   exxfxfx
Método de Newton RaphsonMétodo de Newton-Raphson
E l 8 4 5( ) 0 Exemplo: 8 4.5( sen ) 0x x  
( ) 8 4.5( sen )f x x x   '( ) 4.5(1 cos )f x x  
 Estimativa inicial: x=2
 Erro relativo máximo de 10-3
( ) ( )f ( ) ( )f
 Erro relativo máximo de 10 3
1
8 4.5( sen )
4.5(1 cos )
i i
i i
x xx x
x
    
Método de Newton RaphsonMétodo de Newton-Raphson
i 1 i=1
2
8 4.5(2 sen 2)2 2.48517
4 5(1 cos 2)
x     2 1 2.48517 2 0.24
2
x x  
 i=2
4.5(1 cos 2) 
1 2x
3
8 4.5(2.48517 sen 2.48517)2.48517 2.43099
4.5(1 cos 2.48517)
x     
3 2 2.43099 2.48517 0 02x x 
 i=3
3 2
2
2.43099 2.48517 0.02
2.48517
x x
x
 
4
8 4.5(2.43099 sen 2.43099)2.43099 2.43046
4.5(1 cos 2.43099)
x     
2 43046 2 430993 2
3
2.43046 2.43099 0.0002
2.43099
x x
x
  
Solução: 2.43
Algoritmo do método de NewtonAlgoritmo do método de Newton 
function [res v x it]=fnewton(func dfunc x tol)function [res,v_x,it]=fnewton(func,dfunc,x,tol)
%subrotina para calculo de raiz usando metodo de Newton
%entrada: funcao e derivada analitica
%entrada:x é um valor inicial , tol é a precisao desejada
%saida: res é a raiz v x é o vetor com o valor de x em cada interacao e it%saida: res é a raiz, v_x é o vetor com o valor de x em cada interacao e it
%o numero de interacoes
it=0; x0=x;
v x=[];v_x=[];
d=feval(func,x0)/feval(dfunc,x0);
while abs(d)>tol
v x=[v x x0];
1
( )
'( )
n
n n
n
f xx x
f x
 
v_x=[v_x x0];
x1=x0-d;
it=it+1;
x0=x1;
d=feval(func x0)/feval(dfunc x0);d=feval(func,x0)/feval(dfunc,x0);
end
res=x0;
Algoritmo do método de Newton
function Xs = NewtonRoot(Fun,FunDer,Xest,Err,imax)
% NewtonRoot encontra a raiz de Fun = 0 proxima ao ponto Xest usando o metodo de Newton.
% Variaveis de entrada:
% Fun Nome (string) do arquivo com a funcao que calcula Fun para um dado x.
% FunDir Nome (string) do arquivo com a function que calcula a derivada de
% Fun para um dado x.
% Xest Estimativa inicial da solucao% Xest Estimativa inicial da solucao.
% Err Erro maximo.
% imax Numero maximo e iteracoes
%Variaveis de saida:
% Xs Solucao
for i = 1:imax
Xi = Xest - feval(Fun,Xest)/feval(FunDer,Xest);
if abs((Xi Xest)/Xest) < Errif abs((Xi - Xest)/Xest) < Err
Xs = Xi;
break
end
Xest = Xi;
end
if i == imax
fprintf('Solucao nao foi obtida em %i iteracoes \n' imax)fprintf( Solucao nao foi obtida em %i iteracoes.\n ,imax)
Xs = (‘Sem resposta’);
end
Método da SecanteMétodo da Secante
N t i t ã d i t t Newton: xn+1 interseção do eixo x com a tangente em xn
 Secante: xn+1 interseção do eixo x com a secante da 
curva entre xn-1e xn
xn 1,f(xn 1)
Secante
xn,f(xn)
xn-1,f(xn-1)
raiz
xn-1xn
xn+1
y=f(x)
Método da SecanteMétodo da Secante
( ) ( ) ( ) 0f x f x f x
Inclinação da secante:
1 2 2
1 2 2 3
( ) ( ) ( ) 0f x f x f x
x x x x
  
 1 2
3 2 2( ) ( ) ( )
x x
x x f x
f f
 
Isolando x3:
3 2 2
1 2
( )
( ) ( )
f
f x f x
 1
1
1
( )
( ) ( )
n n
n n n
n n
x x
x x f x
f x f x



  1n n
Método da SecanteMétodo da Secante
A d ti ti i i i d t d l d d As duas estimativas inicias podem estar de um lado da 
solução ou a solução pode estar entre os dois pontos
Método da SecanteMétodo da Secante
 Fórmula recursiva: Fórmula recursiva:
 1
1 ( )
n nx xx x f x 
   1
( ) ( )'( ) n nn
f x f xf x  
 Portanto o método da secante é uma aproximação do
1
1
( )
( ) ( )n n n n n
x x f x
f x f x    1( )n n nf x x 
 Portanto, o método da secante é uma aproximação do 
método de Newton
 Vantagem do método da secante: não é necessário o 
cálculo da derivada analítica
 Desvantagem: necessidade de duas estimativas iniciaises a age ecess dade de duas es a as c a s
Algoritmo do método da SecanteAlgoritmo do método da Secante
function [res v x it]=fsecante1(func x0 x1 tol)function [res,v_x,it]=fsecante1(func,x0,x1,tol)
it=0;
v_x=[];
d=feval(func,x1)*(x1-x0)/(feval(func,x1)-feval(func,x0));
while abs(d)>tol
v x=[v x x1];v_x [v_x x1];
x2=x1-d;
it=it+1;
x0=x1;
1
1
1
( )
( ) ( )
n n
n n n
n n
x xx x f x
f x f x



   
;
x1=x2;
d=feval(func,x1)*(x1-x0)/(feval(func,x1)-feval(func,x0));
end
res=x1; %res=x2 equivalente
it;
Algoritmo do método da Secante
function Xs = SecantRoot(Fun,Xa,Xb,Err,imax)
% SecantRoot encontra a raiz de Fun = 0 usando o metodo da secante.
% Variaveis de entrada:
% Fun Nome (string) do arquivo com a funcao que calcula Fun para dado x% Fun Nome (string) do arquivo com a funcao que calcula Fun para dado x.
% a, b Dois pontos vizinhos a raiz (do mesmo lado da raiz ou em lados opostos).
% Err Erro maximo.
% imax Maximo numero de iteracoes
%Variaveis de saida:
% Xs Solucao
for i = 1:imaxfor i = 1:imax
FunXb = feval(Fun,Xb);
Xi = Xb - FunXb*(Xa-Xb)/(feval(Fun,Xa)-FunXb);
if abs((Xi - Xb)/Xb) < Err(( ) )
Xs = Xi;
break
end
Xa = Xb;Xa = Xb;
Xb = Xi;
end
if i == imax
fprintf('Solucao nao foi obtida em %i iteracoes.\n',imax)
Xs = ('Sem solucao');
end
Exercício 1Exercício 1
E l 3 2 Exemplo:
Z ? N i t õ ?
3 210 29 20 0x x x   
 Zero? No interações?
100
0a) 7x 
50
(
x
)
0)
0b) 2x  
0
f
(
ac) Usar ginput para achar 3 raiz
0 2 4 6 8
0
xx
Exercício 1Exercício 1
El b l it M tl b l l í d li ô iElabore um algoritmo em Matlab que calcule as raízes do polinômio:
3 210 29 20 0x x x   
Pede-se ainda:
 Faça um gráfico da função para mapeá-la. Use o comando ginput para
auxiliar a encontrar a estimativa inicial das raízes.
 Teste diferentes métodos de solução (bisseção, falsa posição, newton,
secante) e compare-os em termos da velocidade de convergência (número
de iterações) e do valor encontrado paraa raiz.
 Teste ainda diferentes estimativas iniciais (por exemplo, 7 e -2),
verificando a sensibilidade do método.
Exercício 1 SoluçãoExercício 1 - Solução
E l 3 2
100
 Exemplo: 3 210 29 20 0x x x   
function F=f302(x);
F=x.^3-10*x.^2+29*x-20; 40
60
80
x
)
f(x)
interacoes de Newton
F x. 3 10 x. 2+29 x 20;
function F=f303(x);
F=3*x.^2-20*x+29; 0
20
40
f
(
x
%calcular raiz com Newton, derivada analítica
[x_N,v_x_N,it_N]=fnewton('f302','f303',7,0.00005);
0 2 4 6 8
-20
x
v_x=0:0.1:8;
plot(v_x,f302(v_x),v_x_N,f302(v_x_N),'ob')
xlabel('x'); ylabel('f(x)');xlabel( x ); ylabel( f(x) );
legend('f(x)','interacoes de Newton')
grid on
fprintf('Raizes de Newton: %6.4f\n',x_N)
it_N
Exercício 1 SoluçãoExercício 1 - Solução
 Exemplo: 3 210 29 20 0x x x   
%calcular raiz com Secante, derivada numérica
[x_sec,v_x_sec,it_sec]=fsecante1('f302',4.5,7,0.00005);
Testar diferentes intervalosfigure
plot(v_x,f302(v_x),v_x_sec,f302(v_x_sec),'ob')
xlabel('x'); ylabel('f(x)'); legend('f(x)','interacoes - Secante')
id
Testar diferentes intervalos
grid on
fprintf('Raizes de secante: %6.4f\n',x_sec)
it sec
100
_
%teste ginput
v_x=0:0.1:8;
40
60
80
(
x
)
f(x)
interacoes - Secante
figure
plot(v_x,f302(v_x))
xlabel('x'); ylabel('f(x)');
legend('f(x)') title('testar ginput')
0
20
40
f
(
legend('f(x)') title('testar ginput')
grid on
aprox_3=ginput(3)
0 1 2 3 4 5 6 7 8
-20
x
Exercício 2Exercício 2
 Repita o que pede o exercício anterior para as seguintes EANL 
  0cos a) 3  xex
(adicionalmente, use o comando fzero do Matlab)
 
  01 c)
01sin b)
5 

x
/x)(
0484410103 d) 2345  xxxxx
 Atenção: explore cada exemplo, imprimindo tabelas com histórico 
de convergência, solução para diferentes estimativas iniciais etc.
Exercício 3
 O fator de fricção de Darcy fD é utilizado em cálculos de
184000
escoamento de fluidos. Para regime turbulento, o fator de fricção
pode ser calculado pela seguinte expressão:
0.5 0.5
Re
1 / 2.51( ) 0.86858ln
3.7D D D
DF f
f N f
     
Re 184000
0.00015 ft
0 17225 ft
N
D



0.17225 ftD 
 onde NRe é o número de Reynolds,  é a rugosidade do da
tubulação e D o diâmetro do tubotubulação e D o diâmetro do tubo.
 Elabore um programa em Matlab que calcule o fator de fricção parap g q ç p
diferentes diâmetros, fazendo um gráfico de fator de fricção versus
diâmetro. No mesmo gráfico, trace curvas para diferentes
rugosidades Atenção para a estimativa inicial!rugosidades. Atenção para a estimativa inicial!
Sistemas não Lineares de Equações AlgébricasSistemas não Lineares de Equações Algébricas
 Qualquer sistema de equações algébricas que não Qualquer sistema de equações algébricas que não 
puder ser expresso na forma matricial típica de um 
sistema linear será, portanto, um SISTEMA NÃO p
LINEAR.
F nção m lti ariá el ( ) 0 1 2f i Função multivariável 1 2( , ,..., ) 0 1,2,...,i nf x x x i n 
 Métodos numéricos (extensão do problema de zero 
de função): Newton-Raphson é o mais utilizado.
 Embora não tenha convergência garantida, é um dos métodos 
mais eficientes, possui uma elevada velocidade de convergência.
Sistemas não Lineares de EASistemas não Lineares de EA
Si t ã li é i    0x,...,x,xf 211 Sistema não linear genérico:   
 



 
0x,...,x,xf
, ,,
n212
n211

 Expansão em série de Taylor
   0x,...,x,xf n21n
     1 1 1 111 1 2 1 1 2 1 1
1
, ,..., 0 , ,...,
k
k k k k k k k k
n n
x x
ff x x x f x x x x x
x
   

     
   1 11 12 2
2
 ...
k k
k k k k
n n
nx x x x
f fx x x x
x x
 
 
        
   1 1 12 1 2 2 1 2, ,..., 0 , ,...,k k k k k kn nf x x x f x x x       12 1 1
1 k
k k
x x
f x x
x
f f


  
    1 12 22 2
2
 ...
k k
k k k k
n n
nx x x x
f fx x x x
x x
 
 
        
Sistemas não Lineares de EASistemas não Lineares de EA
E ã é i d T l Expansão em série de Taylor
       1 1 11 1 11 1 2 1 1 2 20 , ,..., ...k k k k k k k k kn n nf f ff x x x x x x x x xx x x                
       
1 2
1 1 12 2 2
k k knx x x x x x
k k k k k k k k k
x x x
f f ff
  
  
  
         1 1 12 2 22 1 2 1 1 2 2
1 2
0 , ,..., ...
k k k
k k k k k k k k k
n n n
nx x x x x x
f f ff x x x x x x x x x
x x x
  
  
              
 Rearrumando,
       k k k k 1 k k 1 k k 1 k1 1 11 1 2 n 1 1 2 2 n nf f ff x ,x ,...,x x x x x ... x x                     
       
k k k
1 1 2 n 1 1 2 2 n n
1 2 nx x x x x x
k k k k 1 k k 1 k k 1 k2 2 2
2 1 2 n 1 1 2 2 n n
f x ,x ,...,x x x x x ... x x
x x x
f f ff x ,x ,...,x x x x x ... x x
  
  
    
                   
k k k
2 1 2 n 1 1 2 2 n n
1 2 nx x x x x x
, , ,
x x x
 
    
 
Sistemas de EA não linearesSistemas de EA não-lineares
N f t i i l Na forma matricial,
1 1 1f f f     
 
1 2 nk k k k 1
1 1 2 n 1
2 2 2k k k
x x x
f x ,x ,...,x x
f f f
f
k k k   
             
x x x x x x

 
 
k
1
k 1 k
x

    
 
2 2 2k k k
2 1 2 n
1 1 1
k k k
f x ,x ,...,x
x x x
k k k  
                
x x x x x x

    
 
 
k 1 k
2 2
k 1 k
x x     k k kn 1 2 n
n n n
1 1 1
f x ,x ,...,x
f f f
x x x
             x x x x x x

 k 1 kn nx x   
1 1 1k k k   x x x x x x
kF
kJMatriz Jacobiano
1 kk x x
1
1k k k k

   x x J F
kJMatriz Jacobiano
Calcular a cada interação -1kJ
Exercício 4 EquilíbrioExercício 4 - Equilíbrio
P i t bi á i ilíb i LV (lí id ) ã Para um sistema binário em equilíbrio LV (líquido-vapor), a pressão
de cada componente em cada fase é dada por:
0P PLí id 0(1 )P P01 1P xPLíquido:
Vapor: 1P yP 2 (1 )P y P 
0
2 2(1 )P x P 
 onde a pressão de vapor de cada componente na temperatura T é
dada pela equação de Antoine:
1 3848,1A   2 4328,1A  A0 1
1
,
17,532B 
2
2
,
17,913B i
i
i BT
AP 0ln
 No equilíbrio, a pressão da fase vapor é a mesma da fase líquida. A
pressão total de uma fase, por sua vez, é dada pela soma das
pressões parciais:p p
0 0
1 2(1 )P xP x P  
Exercício 4 EquilíbrioExercício 4 - Equilíbrio
P d l b M tl b Pede-se para elaborar um programa em Matlab que:
 Para uma mistura binária equimolar a uma pressão P =760 mmHg,
calcule a temperatura no equilíbrio.p q
 Para a mesma mistura, trace um gráfico T versus P.
 Trace um gráfico contendo curvas T versus composição (x) para
diferentes pressões.
 Para o primeiro item, pelo menos, compare a solução pelos métodos
fzero da secante e Newton (usando derivada analítica e numérica) Afzero, da secante e Newton (usando derivada analítica e numérica). A
derivada analítica da função é dada por:
1 1 2 2
0 0
/ /1 2 1 1 2 2
2 2 2
(1 ) (1 )'( ) A T B A T BxA x A xAP x A Pf T e e
T T T
       
Exercício 5 Equação de Estado (SRK)Exercício 5 – Equação de Estado (SRK)
A ã d t d d S R dli h K (SRK) A equação de estado de Soave-Redlich-Kowng (SRK) para um 
componente pode ser escrita na forma cúbica:
3 2 2( ) 0Z Z A B B Z AB    
 Sendo
( ) 0Z Z A B B Z AB     
2 2
2
T  2 2
0, 4278 0,0867c c
c c
R T R Ta b
P P
  1 1
c
Ts
T
           
2 2 
a P b PA B
R T RT
    48508.055171.115613.0 2  wws
 onde Tc e Pc são a temperatura crítica e a pressão crítica,
respectivamente, w é o fator acêntrico
Exercício 5 Equação de Estado (SRK)Exercício 5 – Equação de Estado (SRK)
El b l it M tl b l l f t d Elabore um algoritmo em Matlab que calcule o fator de
compressibilidade de um componente nas seguintes condições:
500 1 a 40 atm 1 atm = 101325PaT K P 
0,1931
425,2 3799kPa 8314c c
w
T P R

  
 Observe que 3 raízes são encontradas, porém apenas uma é real e
portanto útil.portanto útil.
 De posse do valor de Z, é possível calcular o volume molar ocupado
l ápelo gás:
RT
VPZ


Exercício 5 Equação de Estado (SRK)Exercício 5 – Equação de Estado (SRK)
C f t d ibilid d d ã id lid d d Como o fator de compressibilidade mede a não-idealidade de um
gás, pede-se que o programa:
 Calcule o Z para diferentes temperaturas a uma dada Pp p
 Calcule o volume molar real e ideal
 Trace um gráfico do volume molar real versus ideal para as diferentes Trace um gráfico do volume molar real versus ideal para as diferentes
pressões. No mesmo gráfico, trace uma linha pontilhada de 45º (para
servir de comparação dos pontos).p ç p )
 Trace um gráfico de Z versus T. Adicionalmente, trace uma linha
pontilhada para representar a idealidade (Z=1).
Exercício 6 Equilíbrio QuímicoExercício 6 – Equilíbrio Químico
D õ l l f d C Duas reações ocorrem em paralelo formando C:
1k2A +B C
2kA +D C
 As constantes de equilíbrio de cada reação são dadas por:
45 10CCK  24 10CCK 1 2 5 10C
A B
K
C C
   1 4 10C
A D
K
C C
  
 onde CA, CB, CC e CD são as concentrações de cada reagente no
equilíbrio que podem ser calculadas pelo seguinte balanço de
massa:massa:
Exercício 6 Equilíbrio Químico
B l d
Exercício 6 – Equilíbrio Químico
 Balanço de massa:
,0 1 ,0 2 ,0 1 22 40 30 10
(1 ) 15 15
A A B DC C x C x C x x
C C
     
1 ,0 1
,0 1 ,0 2 ,0 1 2
2 0 2
(1 ) 15 15
15 10
(1 ) 10 10
B B
C C B D
D D
C x C x
C C x C x C x x
C x C x
   
    
   
 Onde x1 e x2 são as conversões de cada reação. As concentrações 
de alimentação de cada reagente são (kmol/m3):
2 ,0 2(1 ) 10 10D DC x C x
de alimentação de cada reagente são (kmol/m ):
,0 ,0 ,0 ,040 15 0 10A B C DC C C C   
 Deseja-se saber quais as conversões de cada reação no equilíbrio.
Exercício 6 Equilíbrio Químico
P t t d l b l it M tl b
Exercício 6 – Equilíbrio Químico
 Para tanto, pede-se para elaborar um algoritmo em Matlab que
calcule as conversões no equilíbrio. Adicionalmente, pede-se:
 Teste diferentes estimativas iniciais e trace um gráfico contendo duasg
curvas: x1 versus estimativa de K1 e x2 versus estimativa de K1. faça o
mesmo para diferentes estimativas de K2. Como se comporta o
gráfico? O que isso significa?
 Mantendo uma estimativa inicial, calcule as conversões para diferentes
concentrações de alimentação de A e trace um gráfico x1 e x2 versus
CA0.
Exercício 7 CombustãoExercício 7 – Combustão
U â d b tã i t d Uma câmara de combustão queima propeno como mostrado na
Figura abaixo.
 Através do balanço de massa atômico é possível calcular a
composição do gás de combustãocomposição do gás de combustão
1 2 P3n n n 
2 3 8n n n
1 2 3 A2 0.21n n n n  
5 A0.79n n
C :
H :
O :
N :
3 4 P2 3 8n n n  5 A0.79n nH : N :
Exercício 7 CombustãoExercício 7 – Combustão
El b l it M tl b l l i ã d á Elabore um algoritmo em Matlab que calcule a composição do gás
de combustão para uma base de cálculo nA = 100 moles/h e
diferentes vazões de ar (2000 a 3000 moles/h). Trace um gráfico da
fração molar de CO e H2O para cada vazão de ar.
 Se a reação for considerada a equaçãoCO + H O H + CO Se a reação for considerada, a equação
para o balanço atômico do N2 deve ser:
CO + H2O  H2 + CO2
 1 4 2 3 0n n K T n n 
 Considerando K =25, resolva o mesmo problema anterior.
 1 4 2 3 0n n n n
 Para uma dada vazão de entrada, compare as duas solução e
comente.