Unicamp - Modelagem - Aula 3

Prévia do material em texto

Aula 3
Sistematização de Estratégias de Solução de 
Sistemas a Equações Algébricas
Porque estudar Sistemas Lineares?
• BME Lineares – etapa importante
– Visão sistêmica – permite o desenvolvimento de algoritmos
– Dimensão do problema – milhares de equações
• Métodos Iterativos – Convergência
• Esparsidade – Métodos especializados
• Estrutura
• Etapa importante na solução de problemas não lineares – métodos de aproximação
– Todos os aspectos anteriores se repetem......
• Simuladores – várias camadas de interpretação das entradas para converter em um 
procedimento numérico....
Lineares
Não Lineares
Tópicos
Lineares: surgem com frequência em problemas de balanço de massa
• Dois em exercícios em sala de aula com objetivo de mostrar situações importantes em que os problemas de EQ são 
modelados como sistemas lineares
• Importante para entender a estrutura de desenvolvimento de um simulador de processos
Não Lineares: surgem em um número expressivo de problemas da EQ
Análise Estrutural
• Ferramenta para a análise de particionamento, variável de partição e ordem de solução de um sistema de equações
Estrategias de solução: exploram a estrutura do problema
Solução simultânea
Como escolher a variável de 
partição? (tearing)
Particionamento
Análise estrutural
O que você entende por solução 
simultânea?
VISÃO SISTÊMICA
Exemplo de motivação: mesmo sendo um 
assunto novo, modelem o balanço de massa da 
coluna de absorção a seguir. Basta escrever 
para cada prato: 
massa entrada=massa saída
Exemplo de aplicação: coluna de absorção com seis pratos
Análise de estrutura
• Deseja-se calcular as composições x6 e y1 à saída de uma 
coluna de absorção. O estudo desta coluna de absorção 
mostrou que pode-se adotar um modelo linear de 
equilíbrio entre as fases do tipo ym=a*xm+b
– Desenvolva o modelo explicitando todas as hipóteses
– Sabendo que a= 0,72; b=0; G=66,7 kmol/min; L=40,8 kmol/min, 
resolva o modelo desenvolvido para as seguintes situações:
• x0=0 e na alimentação de gás se tem 0,2 kmol de soluto/kmol de inerte
• x0=0 e na alimentação de gás se tem 0,3 kmol de soluto/kmol de inerte
Coluna de absorção
O que é uma coluna de absorção?
Dispositivo de contato em 
contracorrente em que há transferência 
de massa de um componente presente 
na fase gasosa para a fase liquida. 
Portanto, a fase liquida será enriquecida 
com este componente
Fase liquida 
enriquecida
Modelo
Quais são as principais 
características que 
aparecem neste problema?
ZEROS Algoritmos especializados
-Estrutura
-Esparsidade
yn=a*xn+b
Sistemas Lineares: Métodos de solução
• Eliminação de Gauss
• Decomposição LU
• Gauss Siedel
Solução de Sistemas
• Dimensão
• Estrutura
• Estratégias que existem para resolver sistemas lineares
– Eliminação de Gauss
– Gauss Siedel (iterativo)
– Decomposição LU
– .....
Eliminação de Gauss – Revisão Rápida
Método de solução de equações [A][X]=[C]
Dois passos
1. Eliminação a frente ( Forward)
2. Retrosubstituição (Backward)
Eliminação










−=




















−−
735.0
21.96
8.106
 
7.000
56.18.40
1525
3
2
1
x
x
x










=




















2.279
2.177
8.106
 
112144
1864
1525
3
2
1
x
x
x
O objetivo da eliminação é transformar a matriz dos coeficientes
em uma matriz triangular superior
Eliminação à frente
1
11
21
1
11
21
212
11
21
121 ... b
a
a
xa
a
a
xa
a
a
xa nn =+++
22323222121 ... bxaxaxaxa nn =++++
1
11
21
21
11
21
2212
11
21
22 ... b
a
a
bxa
a
a
axa
a
a
a nnn −=







−++







−
'
2
'
22
'
22 ... bxaxa nn =++
Eliminação na linha 2
−
_________________________________________________
ou
Algoritmo
Eliminação à frente
No passo (n-1) da eliminação, o Sistema de equações será:
'
2
'
23
'
232
'
22 ... bxaxaxa nn =+++
"
3
"
33
"
33 ... bxaxa nn =++
( ) ( )11 −− =
n
nn
n
nn bxa
. .
. .
. .
11313212111 ... bxaxaxaxa nn =++++
Matriz resultante do procedimento de eliminação à frente
















=
































− )(n-
n
"
'
n
)(n
nn
"
n
"
'
n
''
n
b
b
b
b
x
x
x
x
a
aa
aaa
aaaa
1
3
2
1
3
2
1
1
333
22322
1131211
0000
00
0




Substituição –Passo 2










−=




















−−
735.0
21.96
8.106
 
7.000
56.18.40
1525
3
2
1
x
x
x
Exemplo de um Sistema de 3 equações
Substituição Reversa
'
2
'
23
'
232
'
22 ... bxaxaxa nn =+++
"
3
"
3
"
33 ... bxaxa nn =++
( ) ( )11 −− =
n
nn
n
nn bxa
. .
. .
. .
11313212111 ... bxaxaxaxa nn =++++
Substituição reversa
( ) ( ) ( ) ( )
( ) 1,...,1for
...
1
1
,2
1
2,1
1
1,
1
−=
−−−−
=
−
−
+
−
++
−
+
−
ni
a
xaxaxab
x
i
ii
n
i
nii
i
iii
i
ii
i
i
i
( ) ( )
( ) 1,...,1for1
1
11
−=
−
=
−
+=
−−
ni
a
xab
x
i
ii
n
ij
j
i
ij
i
i
i
)1(
)1(
−
−
=
n
nn
n
n
n
a
b
x
Algoritmo
Exemplo 1
Tempo, Veloc, 
5 106.8
8 177.2
12 279.2
A velocidade de subida de um foguete é medida em três pontos
Os dados podem ser aproximados por uma polinomial:
( ) 12.t5 , 32
2
1 ++= atatatv
Determine a velocidade em t=6 segundos .
( )s t ( )m/s v
Tabela 1 velocidade vs tempo.
Exemplo 1 
( ) 12.t5 ,atatatv ++= 32
2
1
Assumindo










=




















3
2
1
3
2
3
2
2
2
1
2
1
1
1
1
v
v
v
a
a
a
 
tt
tt
tt
3
2
1
A matriz resultante será:
Usando os dados da tabela será










=




















2.279
2.177
8.106
112144
1864
1525
3
2
1
a
a
a
 
Dividir Eq1 por 25 e multiplicá-la por 64
Eliminação à frente










2.279112144
2.1771864
8.1061525



.
   408.27356.28.126456.28.1061525  =
 
 
 208.96 56.18.4 0 
408.273 56.2 8.1264
177.2 1 8 64 
−−−
−













−−−
2.279112144
208.9656.18.40
8.1061525



56.2
25
64
=
Subtraia o resultado da Eq 2
Substitua a nova equação na
equação 2
Eliminação à frente










−−−
2.279112144
208.9656.18.40
8.1061525



.
   168.61576.58.2814476.58.1061525  =
 
 
 968.335 76.48.16 0 
168.615 76.5 8.28 144
279.2 1 12 144 
−−−
−













−−−
−−−
968.33576.48.160
208.9656.18.40
8.1061525



Dividir a equação 1 por 25 e multiplicá-la 
por 144
76.5
25
144
=
Subtrair o resultado da equação 3
Substituir o resultado na eq 3
Eliminação à frente
   728.33646.58.1605.3208.9656.18.40 −−−=−−− 










−−−
−−−
968.33576.48.160
208.9656.18.40
8.1061525



 
 
 .760 7.0 0 0 
728.33646.516.80
335.968 76.416.80 



−−−−
−−










−−−
76.07.000
208.9656.18.40
8.1061525



Dividir a eq 2 por -4,8 e multiplicar
por -16,8 
5.3
8.4
8.16
=
−
−
Subtrair o resultado da eq 3
Substituir a eq 3 pela nova eq 3
Algoritmo – Back Substitution
Exemplo 1- Solução Final
( ) ( ) ( )
.m/s 686.129 
08571.166905.196290472.06
2
=
++=vSolução
O vetor solução é:










=










08571.1
6905.19
290472.0
3
2
1
a
a
a
O polinomial que passa por estes pontos é:
( )
125 ,08571.16905.19290472.0 2
32
2
1
++=
++=
ttt
atatatv
955
11326
3710
321
321
32
=+−
=++
=−
xxx
xxx
xx










=




















−
−
9
11
3
 
515
326
7100
3
2
1
x
x
x
Armadilha – Divisão por zero
Neste caso tem divisão por zero?
955
14356
1571012
321
321
321
=+−
=++
=−+
xxx
xxxxxx










=




















−
−
9
14
15
 
515
356
71012
3
2
1
x
x
x
Tem divisão por zero aqui? SIM
28524
14356
1571012
321
321
321
=+−
=++
=−+
xxx
xxx
xxx










=




















−
−
28
14
15
 
5124
356
71012
3
2
1
x
x
x










−
=




















−
−
2
5.6
15
 
192112
5.600
71012
3
2
1
x
x
x
A divisão por zero é uma possibilidade em qualquer passo
da eliminação de Gauss
Armadilha 2 - Erros de Arredondamento










=










1
1
1
 
3
2
1
x
x
x










=




















−−
9
751.1
45
 
315
7249.23
101520
3
2
1
x
x
x
Solução exata










=










999995.0
05.1
9625.0
 
3
2
1
x
x
x










=










99995.0
5.1
625.0
 
3
2
1
x
x
x
5 algarismos
significativos
6 algarismos
significativos
Exemplo 2










=




















2279
2177
8106
112144
1864
1525
3
2
1
.
.
.
a
a
a
 
Resolva o problema abaixo fazendo o 
pivoteamento parcial
Resolva pelos dois métodos e verifique as respostas










=




















2279
2177
8106
112144
1864
1525
3
2
1
.
.
.
a
a
a
 










=










15.1
67.19
2917.0
3
2
1
a
a
a










=










08571.1
6905.19
290472.0
3
2
1
a
a
a
Com 
Pivoteamento
Parcial
Sem pivoteamento parcial
Método de Gauss-Seidel
Método iterativo (porque usar um método iterativo para resolver um 
Sistema linear??)
Procedimento Básico:
-Resolver algebricamente cada equação com respeito a xi 
-Assumir um vetor de soluções (chute inicial)
-Resolver para cada xi e repetir
-Usar o valor aproximado do modulo do erro relativo para avaliar convergencia.
Método Gauss-Seidel – Porque?
Permite o controle de erros de arredondamento
Se os aspectos físicos da modelagem são bem compreendidos, podem ser feitas boas 
estimativas iniciais, diminuindo o número de iterações necessárias para resolver o problema 
Método de Gauss-Seidel
Algoritmo
Reescrevendo cada equação
11
13132121
1
a
xaxaxac
x nn
−−−
=

nn
nnnnnn
n
nn
nnnnnnnnn
n
nn
a
xaxaxac
x
a
xaxaxaxac
x
a
xaxaxac
x
11,2211
1,1
,122,122,111,11
1
22
23231212
2
−−
−−
−−−−−−−
−
−−−−
=
−−−−
=
−−−
=




Da equação 1
Da equação 2
Da equação n-1
Da equação n
Método de Gauss-Seidel
Algoritmo
Forma geral para qualquer linha ‘i’
.,,2,1,
1
ni
a
xac
x
ii
n
ij
j
jiji
i =
−
=


=
Quando esta equação pode ser usada?
Exemplo 1
Tempo, Veloc, 
5 106.8
8 177.2
12 279.2
A velocidade de subida de um foguete é medida em três pontos
Os dados podem ser aproximados por uma polinomial:
( ) 12.t5 , 32
2
1 ++= atatatv
Determine a velocidade em t=6 segundos usando o Método de Gauss-Siedel
( )s t ( )m/s v
Tabela 1 velocidade vs tempo.
Método de Gauss-Seidel : Exemplo 1










=




















2.279
2.177
8.106
112144
1864
1525
3
2
1
a
a
a
 Estimativa inicial da solução
(Chute inicial) 









=










5
2
1
3
2
1
a
a
a
25
58.106 32
1
aa
a
−−
=
8
642.177 31
2
aa
a
−−
=
1
121442.279 21
3
aa
a
−−
=
6720.3
25
)5()2(58.106
a1 =
−−
=
( ) ( )
8510.7
8
56720.3642.177
a 2 −=
−−
=
( ) ( )
36.155
1
8510.7126720.31442.279
a3 −=
−−−
=
Método de Gauss-Seidel : Exemplo 1










−
−=










36.155
8510.7
6720.3
3
2
1
a
a
a
100
−
=
new
i
old
i
new
i
ia x
xx
%76.72100
6720.3
0000.16720.3
1a
=
−
= x
%47.125100
8510.7
0000.28510.7
2a
=
−
−−
= x
%22.103100
36.155
0000.536.155
3a
=
−
−−
= x
Ao final da primeira iteração
O erro (simplificando a nomenclatura….) 
é 125.47%
Gauss-Seidel Method: Exemplo 1










−
−=










36.155
8510.7
6720.3
3
2
1
a
a
a
Iteração número 2
Usando
( )
056.12
25
36.1558510.758.106
1 =
−−−
=a
( )
882.54
8
36.155056.12642.177
2 −=
−−
=a
( ) ( )
34.798
1
882.5412056.121442.279
3 −=
−−−
=a
Da iteração1
Método de Gauss-Seidel : Exemplo 1
%543.69100
056.12
6720.3056.12
1a
=
−
= x
( )
%695.85100x
882.54
8510.7882.54
2
=
−
−−−
=a
( )
%540.80100
34.798
36.15534.798
3a
=
−
−−−
= x
Ao final da segunda iteração










−
−=










54.798
882.54
056.12
3
2
1
a
a
a
O erro é 85.695%
Iteration a1 a2 a3
1
2
3
4
5
6
3.6720
12.056
47.182
193.33
800.53
3322.6
72.767
69.543
74.447
75.595
75.850
75.906
−7.8510
−54.882
−255.51
−1093.4
−4577.2
−19049
125.47
85.695
78.521
76.632
76.112
75.972
−155.36
−798.34
−3448.9
−14440
−60072
−249580
103.22
80.540
76.852
76.116
75.963
75.931
Método de Gauss-Seidel : Exemplo 1










=










0857.1
690.19
29048.0
a
a
a
3
2
1
%
1a
 %
2a
 %
3a

Os erros não estão decrescendo a taxas importantes
A solução não está convergindo para a solução real
Método de Gauss Siedel: Armadilha
O que deu errado?
Este exemplo muito simples mostra que nem sempre o método de Gauss Siedel converge
Há uma solução?
Há classe de problemas sempre converge: Coeficiente de matriz diagonalmente dominante
Diagonalmente dominante: [A] em [A] [X] = [C] é diagonalmente dominante se:


=

n
j
j
ijaa
i
1
ii 

=

n
ij
j
ijii aa
1
Para todo‘i’ é Para pelo menos um ‘i’ 
O que é importante ao final desta aula?
Não se pretende ser exaustivo, mas alertar para questões importantes:
• Influencia da estrutura
• Influencia das variáveis (natureza e ordem de grandeza)
Buscar sempre softwares de boa qualidade, que usem boas estratégias de 
estabilização numérica de matrizes
• Lembre-se que ao final sempre resolvemos matrizes....
• Não despreze problemas de pequena dimensão
Analisar sempre as soluções encontradas
Balanço Material para a destilação do álcool.
• Uma companhia planeja produzir álcool comercial pelo processo de 
destilação apresentado a seguir. O processo contem duas colunas 
de destilação, ambas com razões de refluxo 3 para 1. A alimentação 
contem 10.000 kg/h com 80% em massa de água, 10% de álcool e 
10% de material orgânico. O destilado da primeira coluna tem 60% 
em álcool enquanto o da segunda coluna possui 95% em peso. A 
corrente de fundo da primeira coluna contem 80% do material 
orgânico alimentado enquanto o restante sai na corrente de fundo 
da segunda coluna. Não há presença de álcool em qualquer das 
correntes de fundo. Calcule, usando uma abordagem global a 
quantidade de cada componente em cada coluna
– Qual a dimensão deste problema?
Balanço de Massa em uma planta de destilação de 
álcool