Método-Newton-Raphson

Prévia do material em texto

Método de Newton Raphson Estendido para
Sistemas não-Lineares
EQ502 – Introdução à Análise de Processos
1 Método de Newton Raphson
Chapra [1] comenta que talvez a fórmula mais amplamente usada para lo-
calizar uma ráız seja a equação de Newton Raphson (ver Figura 1). Se a
aproximação inicial da ráız for xi, pode-se estender uma reta tangente a par-
tir do ponto [xi, f(xi)] (que proporciona a inclinação f
′
(xi) no ponto xi). O
ponto onde essa tangente cruza a abscissa (eixo x) é uma estimativa da ráız.
Aplicando-se sucessivamente esta aproximação, se a estimativa inicial for con-
vergente (normalmente isto é obtido quando a estimativa inicial é próxima
ao valor da ráız), com poucas iterações se obtém uma boa aproximação da
ráız.
O método de Newton Raphson é para uma função com apenas uma
variável independente e é obtido a partir do teorema de Taylor usando-se
apenas o primeiro termo da expansão. Para N termos, o Teorema de Taylor
ao redor de um ponto xi é dado por:
f(xi+∆x) = f(xi)+
(
df
dx
)
i
∆xi+1+
(
d2f
dx2
)
i
∆x2i+1
2!
+
(
d3f
dx3
)
i
∆x3i+1
3!
+. . . (1)
Nesse caso, ∆xi+1 = xi+1 − xi
Simplificando
f(xi+1) = f(xi)+
(
df
dx
)
i
(xi+1−xi)+
(
d2f
dx2
)
i
(xi+1 − xi)2
2!
+
(
d3f
dx3
)
i
(xi+1 − xi)3
3!
+. . .
(2)
Utilizando apenas o primeiro termo da expansão:
f(xi+1) = f(xi) +
(
df
dx
)
i
(xi+1 − xi) (3)
1
Figura 1: Representação gráfica do método de Newton Raphson (essa figura
é do livro Chapra-Canale)
Como se deseja a ráız da função f(x), faz-se que f(xi+1) = 0. Obtém-se
então:
0 = f(xi)+
(
df
dx
)
i
(xi+1 − xi) ou
(
df
dx
)
i
(xi+1 − xi) = −f(xi) (4)
xi+1 = xi −
f(xi)(
df
dx
)
i
(5)
Notem novamente que o valor em i + 1 foi obtido somente em função dos
valores já conhecidos em i. A Figura 1 mostra graficamente o uso dessa
equação após uma aproximação.
2 Propriedades de matrizes
O método de Newton Raphson Estendido leva a um sistema matricial. Sendo
assim, é interessante relembrarmos algumas propriedades de matrizes.
2
Matriz identidade A matriz identidade é definida por
I =

1 0 . . . 0
0 1 . . . 0
...
...
. . . 0
0 0 . . . 1
 (6)
Dado uma matriz [A], podemos escrever que
[A] [I] = [I] [A] = [A] (7)
Matriz inversa O inverso da matriz A é representado por [A]−1 respeitando
a seguinte restrição:
[A] [A]−1 = [A]−1 [A] = [I] (8)
Resolução de um sistema matricial Dado o seguinte sistema matricial
[A] [X] = [B] (9)
Podemos multiplicar ambos os lados da Equação 9 por [A]−1
[A] [A]−1 [X] = [A]−1 [B] (10)
Chegamos que
[X] = [A]−1 [B] (11)
3 Método de Newton Raphson Estendido
O Método de Newton Raphson Estendido é uma aproximação para a ob-
tenção de uma ráız para uma função de inúmeras variáveis independentes.
Este método é particularmente importante na determinação de parâmetros
e/ou variáveis que correspondem à solução numérica de um conjunto de
Equações que formam um sistema não-linear. Suponhamos que tenhamos
um sistema com três tanques interligados, assim como mostra a Figura 2.
Neste caso não há reação qúımica e deseja-se saber o valor das concentrações
CT1 , CT2 e CT3 . Neste caso estamos considerando apenas três variáveis, mas
o procedimento é facilmente estendido para N variáveis. O balanço de massa
em regime permanente para os três tanques fornece as seguintes equações:
3
Figura 2: Sistema com três tanques interligados
Tanque 1 −→ Q1 × CA01 − Q12 × CT1 + Q31 × CT3 = 0
Tanque 2 −→ Q12 × CT1 + Q2 × CA02 − Q23 × CT2 = 0
Tanque 3 −→ Q23 × CT2 − Q31 × CT3 + Q33 × CT3 = 0
(12)
Substituindo os valores numéricos
Tanque 1 −→ 5× 10 − 7× CT1 + 2× CT3 = 0
Tanque 2 −→ 7× CT1 + 1× 1 − 8× CT2 = 0
Tanque 3 −→ 8× CT2 − 2× CT3 + 6× CT3 = 0
(13)
Rearranjando
Tanque 1 −→ −7× CT1 + 2× CT3 = −50
Tanque 2 −→ 7× CT1 − 8× CT2 = −1
Tanque 3 −→ 8× CT2 − 8× CT3 = 0
(14)
As Equações 14 são solucionadas por se resolver o sistema linear abaixo: −7 0 27 −8 0
0 8 −8
 ×
 CT1CT2
CT3
 =
 −50−1
0
 (15)
Chamando  a1,1 a1,2 a1,3a2,1 a2,2 a2,3
a3,1 a3,2 a3,3
 =
 −7 0 27 −8 0
0 8 −8
 (16)
4
e  a1,4a2,4
a3,4
 =
 −50−1
0
 (17)
O sistema linear acima pode ser solucionado utilizando-se o teorema de
Gauss. Chamemos de
Linha 1 =
Linha 2 =
Linha 3 =
 a1,1 a1,2 a1,3a2,1 a2,2 a2,3
a3,1 a3,2 a3,3
 ×
 xy
z
 =
 a1,4a2,4
a3,4
 (18)
Alternativamente
Linha 1 −→ a1,1 x + a1,2 y + a1,3 z = a1,4 (19)
Linha 2 −→ a2,1 x + a2,2 y + a2,3 z = a2,4 (20)
Linha 3 −→ a3,1 x + a3,2 y + a3,3 z = a3,4 (21)
Usa-se a Equação da Linha 1 para eliminar x das Equações das Linhas
2 e 3. Para tal, o primeiro procedimento é multiplicar a Equação Linha 1
por a2,1/a1,1 e subtrair da Equação Linha 2. Analogamente, multiplique a
Equação Linha 1 por a3,1/a1,1 e subtraia da Equação Linha 3. Desta forma
as Equações linhas 2 e 3 se tornam
Linha 1 −→ a1,1 x + a1,2 y + a1,3 z = a1,4 (22)
Linha 2 −→
[
−a2,1 +
a2,1
a1,1
a1,1
]
x +
[
−a2,2 +
a2,1
a1,1
a1,2
]
y+[
−a2,3 +
a2,1
a1,1
a1,3
]
z =
[
−a2,4 +
a2,1
a1,1
a1,4
]
(23)
Linha 3 −→
[
−a3,1 +
a3,1
a1,1
a1,1
]
x +
[
−a3,2 +
a3,1
a1,1
a1,2
]
y+[
−a3,3 +
a3,1
a1,1
a1,3
]
z =
[
−a3,4 +
a3,1
a1,1
a1,4
]
(24)
Desta forma, as linhas 1, 2 e 3 ficam
Linha 1 −→ a1,1 x + a1,2 y + a1,3 z = a1,4 (25)
5
Linha 2 −→
[
−a2,2 +
a2,1
a1,1
a1,2
]
y+[
−a2,3 +
a2,1
a1,1
a1,3
]
z =
[
−a2,4 +
a2,1
a1,1
a1,4
]
(26)
Linha 3 −→
[
−a3,2 +
a3,1
a1,1
a1,2
]
y+[
−a3,3 +
a3,1
a1,1
a1,3
]
z =
[
−a3,4 +
a3,1
a1,1
a1,4
]
(27)
Chamando
a
′
2,j = −a2,j +
a2,1
a1,1
a1,j j = 2, 3, 4 (28)
a
′
3,j = −a3,j +
a3,1
a1,1
a1,j j = 2, 3, 4 (29)
O sistema agora toma a forma
Linha 1 =
Linha 2
′
=
Linha 3
′
=
 a1,1 a1,2 a1,3a′2,2 a′2,3
a
′
3,2 a
′
3,3
 ×
 xy
z
 =
 a1,4a′2,4
a
′
3,4
 (30)
A linha 3
′
é novamente modificada para eliminar a
′
3,2. Usa-se a Equação
da linha 2
′
para eliminar y da Equação da Linha 3. Para tal, o primeiro
procedimento agora é multiplicar a Equação Linha 2
′
por a
′
3,2/a
′
2,2 e subtrair
da Equação Linha
′
3. Desta forma a Equação linha
′
3 se torna.
[
−a′3,2 +
a
′
3,2
a
′
2,2
a
′
2,2
]
y +
[
−a′3,3 +
a
′
3,2
a
′
2,2
a
′
3,3
]
z =
[
−a′3,4 +
a
′
3,2
a
′
2,2
a
′
3,4
]
(31)
Ou
Linha
′
3 −→
[
−a′3,3 +
a
′
3,2
a
′
2,2
a
′
3,3
]
z =
[
−a′3,4 +
a
′
3,2
a
′
2,2
a
′
3,4
]
(32)
Chamando
a
′′
3,j = −a
′
3,j +
a
′
3,2
a
′
2,2
a3,j j = 3, 4 (33)
Chega-se ao sistema
6
Linha 1 =
Linha 2
′
=
Linha 3
′′
=
 a1,1 a1,2 a1,3a′2,2 a′2,3
a
′′
3,3
 ×
 xy
z
 =
 a1,4a′2,4
a
′′
3,4
 (34)
Da Linha 3
′′
obtém-se o valor de z
z =
a
′′
3,4
a
′′
3,3
(35)
Com z obtém-se o valor de y na linha
′
2
y =
a
′
2,4 − z × a
′
2,3
a
′
2,2
(36)
Finalmente x é obtido pela Linha 1
x =
a1,4 − y × a1,2 − z × a1,3
a1,1
(37)
Para o exemplo do tanque veja se chega aos mesmos valores conforme
mostrado abaixo. Após a primeira substituição, para eliminar CT1 partindo
do sistema de Equações 15 −7 0 20 8 −2
0 −8 8
 ×
 CT1CT2
CT3
 =
 −5051
0
 (38)
Após a segunda substituição: −7 0 20 8 −2
0 0 −6
 ×
 CT1CT2
CT3
 =
 −5051
−51
 (39)
Da terceira linha da matriz obtém-se o valor de CT3
CT3 =
−51
−6
= 8, 5 (40)
Com z obtém-se o valor de CT2 na linha
′
2
CT2 =
51 − CT3 × −2
8
=
51 − 8, 5 × −2
8
= 8, 5 (41)
7
Finalmente CT1 é obtido pela Linha 1
CT1 =
−50 − 8, 5 × 0 − 8, 5 × 2
−7
= 9, 57 (42)
O sistema linear mostrado para apenas três equações é estendido para N
equações usando-se o mesmo procedimento. Constroe-se a matriz diagonal e
inicia-se a determinação das variáveis da Nésima Equação até a primeira.
Conforme já mencionado, o Método de Newton Raphson Estendido é uma
aproximação para a obtenção de uma ráız para uma função de inúmeras
variáveis independentes quando o sistema é não-linear. Suponhamos asEquações 43 - 45 referentes aos exerćıcios 6.7 e 6.8 do livro do Fogler para a
obtenção do m-xileno a partir do mesitileno, que será mostrada no Exerćıcio
4.
CH − 0, 021 = −
(
k1C
1/2
H CM + k2C
1/2
H CX
)
τ (43)
CM − 0, 0105 = − k1C1/2H CM τ (44)
CX = −
(
k1C
1/2
H CM − k2C
1/2
H CX
)
τ (45)
O sistema acima não é linear e não pode ser resolvido sem que o mesmo
não seja primeiramente linearizado. A forma de linearização deste sistema é
aproximar as Equações 43 - 45 com o uso do Teorema de Taylor. Mostrando
para um sistema de três equações, suponhamos que se deseje conhecer as
ráızes x, y e z das funções f1(x, y, z) = g1(x, y, z), f2(x, y, z) = g2(x, y, z) e
f3(x, y, z) = g3(x, y, z). Primeiramente constroem-se as restrições
F1(x, y, z) = f1(x, y, z) − g1(x, y, z) = 0 (46)
F2(x, y, z) = f2(x, y, z) − g2(x, y, z) = 0 (47)
F3(x, y, z) = f3(x, y, z) − g3(x, y, z) = 0 (48)
As funções restrições sempre tem que ser igual a zero. Crie as funções
restrições das Equações 43 - 45 e compare seu resultado com as Equações 83
- 85
Expandem-se então as Equações 46 - 48 ao redor de xi, yi e zi com apenas
o primeiro termo da expansão. Desta forma
τ1(xi + ∆xi+1, yi + ∆yi+1, zi + ∆zi+1) = F1(xi, yi, zi) +(
∂F1(x, y, z)
∂x
)
i
∆xi+1 +
(
∂F1(x, y, z)
∂y
)
i
∆yi+1 +
(
∂F1(xi, yi, z)
∂z
)
i
∆zi+1 (49)
8
τ2(xi + ∆xi+1, yi + ∆yi+1, zi + ∆zi+1) = F2(xi, yi, zi) +(
∂F2(x, y, z)
∂x
)
i
∆xi+1 +
(
∂F2(x, y, z)
∂y
)
i
∆yi+1 +
(
∂F2(x, y, z)
∂z
)
i
∆zi+1 (50)
τ3(xi + ∆xi+1, yi + ∆yi+1, zi + ∆zi+1) = F3(xi, yi, zi) +(
∂F3(x, y, z)
∂x
)
i
∆xi+1 +
(
∂F3(x, y, z)
∂y
)
i
∆yi+1 +
(
∂F3(x, y, z)
∂z
)
i
∆zi+1 (51)
onde
∆xi+1 = xi+1 − xi (52)
∆yi+1 = yi+1 − yi (53)
∆zi+1 = zi+1 − zi (54)
Chamando
(F1x)i =
(
∂F1(x, y, z)
∂x
)
i
; (F1y)i =
(
∂F1(x, y, z)
∂y
)
i
; (F1z)i =
(
∂F1(x, y, z)
∂z
)
i
(55)
(F2x)i =
(
∂F2(x, y, z)
∂x
)
i
; (F2y)i =
(
∂F2(x, y, z)
∂y
)
i
; (F2z)i =
(
∂F2(x, y, z)
∂z
)
i
(56)
(F3x)i =
(
∂F3(x, y, z)
∂x
)
i
; (F3y)i =
(
∂F3(x, y, z)
∂y
)
i
; (F3z)i =
(
∂F3(x, y, z)
∂z
)
i
(57)
Substituindo
τ1(xi+1, yi+1, zi+1) = F1(xi, yi, zi))+ (F1x)i (xi+1−xi)+ (F1y)i (yi+1− yi)+ (F1z)i (zi+1− zi)
(58)
τ2(xi+1, yi+1, zi+1) = F2(xi, yi, zi)+ (F2x)i (xi+1−xi)+ (F2y)i (yi+1− yi)+ (F2z)i (zi+1− zi)
(59)
τ3(xi+1, yi+1, zi+1) = F3(xi, yi, zi)+ (F3x))i (xi+1−xi)+ (F3y)i (yi+1− yi)+ (F3z)i (zi+1− zi)
(60)
Da mesma forma que para o método de Newton Raphson, faz-se que
τ1(xi+1, yi+1, zi+1) = 0 (61)
τ2(xi+1, yi+1, zi+1) = 0 (62)
τ3(xi+1, yi+1, zi+1) = 0 (63)
para que as ráızes múltiplas correspondentes às determinações de xi+1, yi+1
e zi+1 possam ser obtidas. Desta forma
9
F1(xi, yi, zi)) + (F1x)i (xi+1 − xi) + (F1y)i (yi+1 − yi) + (F1z)i (zi+1 − zi) = 0
(64)
F2(xi, yi, zi) + (F2x)i (xi+1 − xi) + (F2y)i (yi+1 − yi) + (F2z)i (zi+1 − zi) = 0
(65)
F3(xi, yi, zi) + (F3x))i (xi+1 − xi) + (F3y)i (yi+1 − yi) + (F3z)i (zi+1 − zi) = 0
(66)
Colocando-se de uma forma matricial
 (F1x)i (F1y)i (F1z)i(F2x)i (F2y)i (F2z)i
(F3x)i (F3y)i (F3z)i

i
×
 (xi+1 − xi)(yi+1 − yi)
(zi+1 − zi)

i+1
= −
 F1(x, y, z)F2(x, y, z)
F3(x, y, z)

i
(67)
Os novos valores de xi+1, yi+1 e zi+1 são iterados até que o método convirja.
Para o caso onde xi+1 ̸= 0, yi+1 ̸= 0 e zi+1 ̸= 0
(xi+1 − xi)
xi+1
< ε (68)
(yi+1 − yi)
yi+1
< ε (69)
(zi+1 − zi)
zi+1
< ε (70)
Se xN+1 = 0, yN+1 = 0, zN+1 = 0
∆xi+1 = xi+1 − xi ≤ ϵ (71)
∆yi+1 = yi+1 − yi ≤ ϵ (72)
∆zi+1 = zi+1 − zi ≤ ϵ (73)
Se xi+1 ̸= 0, yi+1 ̸= 0 e zi+1 ̸= 0, um bom valor para a convergência
é ε = 10−3. Se xN+1 = 0, yN+1 = 0 e zN+1 = 0, um bom valor para
a convergência é ϵ = 10−6. O método deve ser iterado até que todas as
variáveis tenham alcançado a convergência. Para a próxima iteração basta
lembrar que
∆xi+1 = xi+1 − xi ou xi+1 = xi +∆xi+1 (74)
∆yi+1 = yi+1 − yi ou yi+1 = yi +∆yi+1 (75)
∆zi+1 = zi+1 − zi ou zi+1 = zi +∆zi+1 (76)
10
3.1 Exerćıcio 4
Hidrodesalquilação de mesitileno em um CSTR
(exerćıcios 6.7 e 6.8 do livro do Fogler, 3ª Edição)
A produção de m-xileno por hidrodesalquilação de mesitileno sobre um ca-
talizador Houndry Detrol 1 envolve as seguintes reações:
Mesitileno (M) + Hidrogênio (H) −→ m-Xileno (X) + Metano (Me)
O m-Xileno também pode sofrer hidrodesalquilação para formar tolueno:
m-Xileno (X) + Hidrogênio (H) −→ Tolueno (T) + Metano (Me)
Onde os sub́ındices são:
M= mesitileno, X= m-xileno, T=tolueno, Me= metano e H= hidrogênio.
A segunda reação é indesejada, porque o m-Xileno é vendido por um
preço mais alto do que o tolueno. Vemos, portanto, que existe uma grande
motivação para se maximizar a produção de m-Xileno.
A hidrogenação do Mesitileno deve ser conduzida isotermicamente a 1500
°R e 35 atm em um reator CSTR na qual a alimentação é de 66,7 mol% de
hidrogênio e 33,3 mol% de mesitileno. As leis de velocidade para as reações
1 e 2 são, respectivamente,
− r1M = k1CM C
1/2
H
r2T = k2CX C
1/2
H
A 1500°R as velocidade espećıficas de reação são:
Reação 1: k1 = 55, 20 (ft
3/(lbmol)0,5/h Reação 2: k2 = 30, 20 (ft
3/(lbmol)0,5/h
Calcule a conversão de hidrogênio e mesitileno, juntamente com as concen-
trações de sáıda de mesitileno, hidrogênio e xileno na sáıda do CSTR. Resolva
o problema usando um tempo espacial τ = V
υ0
= 238ft
3
476ft3/h
= 0, 5h
1Ind. Eng. Chem. Process Dev., 4, 92 (1965), 146 (1966)
11
As concentrações iniciais são:
CH0 =
yH0 P0
RT
=
0, 667 35
0, 73 1500
= 0, 021lbmol/ft3
CM0 =
1
2
CH0 = 0, 0105
e
CX0 = 0
Resolva com o uso do Método de Newton Raphson Estendido.
Façam o balanço de massa. As equações finais que vocês devem chegar
são:
CH − CHo = −
(
k1C
1/2
H CM + k2C
1/2
H CX
)
τ (77)
CM − CM0 = − k1C
1/2
H CM τ (78)
CX =
(
k1C
1/2
H CM − k2C
1/2
H CX
)
τ (79)
Substituindo chega-se nas equações 43 – 45 repetidas abaixo por convêniência:
CH − 0, 021 = −
(
k1C
1/2
H CM + k2C
1/2
H CX
)
τ (80)
CM − 0, 0105 = − k1C1/2H CM τ (81)
CX =
(
k1C
1/2
H CM − k2C
1/2
H CX
)
τ (82)
Construindo-se as funções restrição
F1 = CH − 0, 021 +
(
k1C
1/2
H CM + k2C
1/2
H CX
)
τ = 0 (83)
F2 = CM − 0, 0105 + k1C1/2H CM τ = 0 (84)
F3 = CX +
(
k2C
1/2
H CX − k1C
1/2
H CM
)
τ = 0 (85)
A resolução do problema acima fica como exerćıcio com o uso da ferra-
menta Excell. A entrega do exerćıcio é para o dia 07 de maio até as 18:30 h.
No dia 07 de maio de vem ser entregues apenas dois arquivos:
1 - arquivo único em word ou pdf contendo a resolução do problema.
2 - planilha em Excell utilizada para resolver o problema.
12
Referências
[1] Chapra, S. C. e Canale, R. P., Métodos Numéricos para Engenharia
5a Edição, McGraw Hill, 2008.
[2] Jenson, V. G. and Jeffreys, G. V., Mathematical Methods in Chemical
Engineering, 2nd Edition, Academic Press, 1977.
[3] Davis. M. E., Numerical Methods and Modeling for Chemical Engi-
neers, .
[4] Ruggiero, M. A. and Rocha Lopes, V. L. Cálculo Numérico – Aspec-
tos Teóricos e Computacionais, McGrawhill, 1988.
[5] Boyce, W. E. and DiPrima, R. C. Elementary Differential Equa-
tions and Boundary Value Problems, John Wiley & Sons, Inc,
4thEdition,1992.
13