APOSTILA_ALUNO

Prévia do material em texto

1 
 
UNIVERSIDADE FEDERAL RURAL DO RIO DE JANEIRO 
INSTITUTO DE TECNOLOGIA 
DEPARTAMENTO DE ENGENHARIA QUÍMICA 
 
 
 
 
 
IT 383 – MÉTODOS COMPUTACIONAIS APLICADOS I 
PROFESSORA: MÁRCIA PEIXOTO VEGA DOMICIANO 
 
 
 
 
AULAS 
1. Aula 1 – Programação aplicada à resolução de problemas da Engenharia Química 
 
Queridos alunos e alunas, 
Tudo bem? 
Vamos usar em nossas aulas o programa computacional livre Octave, que pode ser baixado 
através do site: https://www.gnu.org/software/octave/index. 
Conforme reportado em https://pt.wikipedia.org/wiki/GNU_Octave#cite_note-1, “GNU 
Octave é uma linguagem computacional, desenvolvida para computação matemática. Possui 
uma interface em linha de comando para a solução de problemas numéricos, lineares e não-
lineares, também é usada em experimentos numéricos. Faz parte do projeto GNU, é um 
software livre sob os termos da licença GPL. Foi escrito por John W. Eaton. Possui 
compatibilidade com MATLAB, possuindo um grande número de funções semelhantes”. 
Por favor, para a primeira aula, todos com o programa instalado em seus dispositivos 
eletrônicos. 
Nesta primeira aula introdutória, vamos fazer uma revisão sobre matrizes. 
Revisão sobre matrizes 
Determinante: 
Vamos calcular o determinante de uma matriz A (operação definida somente para 
matrizes quadradas): 













44434241
34333231
24232221
14131211
aaaa
aaaa
aaaa
aaaa
A 
O determinante da matriz resultante da eliminação da linha 2 (i) e da coluna 1(j) da 
matriz A é chamado 21M ou ijM : 
3 
 
444342
343332
141312
ij21
aaa
aaa
aaa
MM  
Vamos definir o cofator ijA : 
  ijjiij M1A  
Vamos efetuar a expansão de Laplace para o cálculo do determinante: 



n
1k
ikikAaA ou 


n
1k
kjkjAaA 
Exemplo de uma matriz 2 por 2: 







2221
1211
aa
aa
A 
Vamos fazer a expansão da linha 1: 


2
1k
k1k1 AaA 
Vamos calcular os cofatores: 
  22221111 aa1A   
  21212112 aa1A   
Logo, 2112221112121111 aaaaAaAaA  que é a velha fórmula de multiplicar 
em cruz. 
Regra de Cramer: 
Vamos relembrar a regra de Cramer: 
Para o sistema: bAx  , a solução é: 
 
 Adet
Adet
x
j
j  , sendo: 



















nn1j,nj,n
n21j,2j,2
n11j,1j,1
1j,n2n1n
1j22221
1j,11211
aaa
aaa
aaa
aaa
aaa
aaa
A








;



















nn1j,nn
n21j,22
n11j,11
1j,n421n
1j22221
1j,11211
j
aab
aab
aab
aaa
aaa
aaa
A








. 
1.1. Comandos 
As linhas de programa começadas por % constituem comentários ou observações. 
Uma linha começando com ‘%‘ é utilizada para armazenar os comentários e as observações 
do programador, e tais comentários e observações não são executados. Se os comentários e 
observações necessitarem de mais de uma linha de programa, cada linha deve começar com 
%. 
O ponto-e-vírgula é utilizado para suprimir a impressão na tela do monitor. Se o 
último caractere de uma sentença for um ponto-e-vírgula, a impressão é suspensa; o 
comando, contudo, é executado, mas o resultado não é mostrado. Além disso, quando se 
entra com os elementos de uma matriz, utiliza-se um ponto-e-vírgula para indicar o fim de 
cada linha, exceto a última. 
/ É o operador divisão direita. B/A é uma divisão escalar ou matricial (equivalente a 
B*inv(A)). 
+ É o operador mais. X + Y efetua a adição das matrizes X e Y, que devem apresentar a 
mesma dimensão, a não ser que uma das variáveis seja um escalar. 
- É o operador menos. X - Y efetua a subtração das matrizes X e Y. X e Y devem 
apresentar a mesma dimensão, a não ser que uma das variáveis seja um escalar. 
5 
 
* É o operador multiplicação matricial. X*Y é o produto entre as matrizes X e Y. Um 
escalar (matriz de 1 por 1) pode ser multiplicado por outro escalar ou por uma matriz. Para o 
caso de multiplicação de matrizes, o número de colunas de X deve ser igual ao número de 
linhas de Y. 
^ É o operador potência. Z = X^y calcula X elevado à potência de y. 
O comando det(A) fornece o determinante da matriz quadrada A. 
Uma característica conveniente do OCTAVE é que as variáveis não precisam ser 
dimensionadas antes do uso. Todos os cálculos no OCTAVE são executados em precisão 
dupla. Contudo, a visualização pode ser com vírgula fixa (ponto fixo) e quatro casas decimais 
- format short. Se pelo menos um dos elementos de uma matriz não for um inteiro exato, há 
quatro formatos de saída possíveis. A saída mostrada na tela pode ser controlada por meio 
dos seguintes comandos principais: 
format short (formato curto) 
format long (formato longo) 
format short e (formato curto com notação científica) 
format long e (formato longo com notação científica) 
Sempre que o OCTAVE for acionado sem nenhum comando específico de formato de 
visualização, o OCTAVE apresentará o formato short. 
No OCTAVE, as variáveis são geradas automaticamente ao serem utilizadas. (As 
dimensões das variáveis podem ser alteradas mais tarde, se necessário.) Tais variáveis 
permanecem na memória até que se entre com um dos comandos exit ou quit. Para se obter 
uma lista das variáveis na área de trabalho, digita-se simplesmente o comando who. Em 
seguida, todas as variáveis correntemente na área de trabalho aparecem na tela do monitor. 
O comando clear elimina da área de trabalho todas as variáveis não-permanentes. 
Para limpar somente uma variável particular, por exemplo ‘x’, o comando clear x deve ser 
empregado. 
Vamos escrever um programa em OCTAVE seguindo as etapas: 
Clicar duas vezes com o mouse no ícone de programa OCTAVE 
Clicar em FILE;NEW;M-FILE 
Aparece uma janela tipo notepad 
Vamos salvar o arquivo: 
Clicar em ARQUIVO; SALVAR COMO 
Onde aparece SALVAR EM vamos procurar o diretório que desejamos salvar o 
arquivo. 
Onde aparece NOME DO ARQUIVO vamos digitar rpeq1.m 
Onde aparece SALVAR COMO TIPO, vamos selecionar TODOS OS ARQUIVOS 
(*.*). 
Para rodar um programa OCTAVE, devemos digitar o nome do programa sem a 
extensão na janela de comando do OCTAVE e apertar o ENTER no teclado. 
 
 
 
 
 
 
7 
 
1.2. Exemplo de programas no OCTAVE 
1.2.1. Programa rpeq1.m 
Um gás natural apresenta a seguinte composição em % molar: 
Componente i % molar ( ix ) 
Peso molecular ( iPM , 
kgmol/kg ) 
4CH 75% 16 
62HC 20% 30 
2N 5% 28 
 
Adotando-se uma base de 100 kgmol, a massa do componente i é: iii PMnm  , 
sendo: in número de moles do componente i. Portanto, a composição ponderal será: 


i
i
i
i m
m
100ponderal% . 
Vamos calcular a composição ponderal e o peso molecular médio, sendo 
 
i
ii xPMPM . 
Início do programa a ser digitado: 
% Um gás natural apresenta a seguinte composição em % molar: 
% Componente 1:CH4 - 75% (PM=16) 
% Componente 2:C2H6 - 20% (PM=30) 
% Componente 3:N2 - 5% (PM=28) 
% Calcular a composição ponderal e o peso molecular médio 
% Peso molecular dos componentes 
PM1=16; 
PM2=30; 
PM3=28; 
% Composição molar dos componentes 
X1=0.75; 
X2=0.20; 
X3=0.05; 
% Cálculo da composição ponderal adotando-se uma base de 100 kgmol 
% O número de moles dos componentes 1, 2 e 3 será: 
n1=100*0.75; 
n2=100*0.20; 
n3=100*0.05; 
% A massa de cada componente será: 
m1=n1*PM1; 
m2=n2*PM2; 
m3=n3*PM3; 
% A massa total será: 
MT=m1+m2+m3; 
% A composição ponderal será: 
pcentpond1=100*(m1/MT) 
pcentpond2=100*(m2/MT) 
pcentpond3=100*(m3/MT) 
% Cálculo do peso molecular médio 
PMMEDIO=PM1*X1+PM2*X2+PM3*X3 
 
Fim do programa a ser digitado 
 
1.2.2. Programa comando1.m 
Início do programa a ser digitado: 
 
%(formato curto - Scaled fixed point format with 5 digits.) 
format short 
X1=[0 1 2 3 4 5 6 7 8 9 10] 
Y1=rand(1,11) 
X2=[0 1 2 3 4 5] 
Y2=rand(1,6) 
pause 
%(formato longo - Scaled fixed point format with 15 digits.) 
format long 
X1=[0 1 2 3 4 5 6 7 8 9 10] 
Y1=rand(1,11) 
X2=[0 1 2 3 4 5] 
Y2=rand(1,6) 
pause 
%(formato curto comnotação científica - 
%Floating point format with 5 digits.) 
format short e 
X1=[0 1 2 3 4 5 6 7 8 9 10] 
Y1=rand(1,11) 
X2=[0 1 2 3 4 5] 
9 
 
Y2=rand(1,6) 
pause 
%(formato longo com notação científica - 
%Floating point format with 15 digits.) 
format long e 
X1=[0 1 2 3 4 5 6 7 8 9 10] 
Y1=rand(1,11) 
X2=[0 1 2 3 4 5] 
Y2=rand(1,6) 
format short 
pause 
who 
pause 
clear X1 
pause 
who 
pause 
X1=[0 1 2 3 4 5 6 7 8 9 10]; 
Y1=rand(1,11); 
X2=[0 1 2 3 4 5]; 
Y2=rand(1,6); 
pause 
who 
pause 
Z=[1 2 3;4 5 6;7 8 9]; 
pause 
Z=[1 2 3;4 5 6;7 8 9] 
pause 
Z=[rand(10,7)] 
pause 
Z(1,:) 
pause 
Z(:,1) 
pause 
X= 1.234 + 2.345 + 3.456 + 4.567 + 5.678 + 6.789... 
+ 7.890 + 8.901 - 9.012 
X= 1.234 + 2.345 + 3.456 + 4.567 + 5.678 + 6.789..... 
+ 7.890 + 8.901 - 9.012 
 
 
 
 
 
 
 
 
1.3. Exercício proposto: Programa rpeq3.m 
Programa repq3.m 
Um gás natural apresenta a seguinte composição em % molar: 
 
Componente i % molar ( ix ) 
Peso molecular ( iPM , 
kgmol/kg ) 
4CH 75% 16 
62HC 20% 30 
2N 2.5% 28 
83HC 2.5% 44 
 
Adotando-se uma base de 100 kgmol, a massa do componente i é: iii PMnm  , 
sendo: in número de moles do componente i. Portanto, a composição ponderal será: 


i
i
i
i m
m
100ponderal% . 
Calcular a composição ponderal e o peso molecular médio, sendo  
i
ii xPMPM . 
11 
 
2. Aula 2 – Programação aplicada à resolução de problemas da Engenharia Química 
Queridas alunas e alunos, 
Queridos alunos e alunas, 
Tudo bem? 
Vamos resolver um sistema de equações algébricas lineares usando a regra de Cramer e 
aprender a traçar gráficos 2D e 3D. 
2.1. Comandos 
O operador dois-pontos pode ser utilizado para criar vetores, indexar matrizes e para 
especificar iterações for. Por exemplo, i:k é o mesmo que [i i + 1 ... k], A(:, j) é a j-ésima 
coluna da matriz A e A (i, :) é a i-ésima linha de A. 
Se a sentença a ser digitada for muito longa, com comprimento superior ao de uma 
linha, podem ser utilizadas reticências constituídas por três ou mais pontos(...), um exemplo 
é: 
x = 1.234 + 2.345 + 3.456 + 4.567 + 5.678 + 6.789... 
+ 7.890 + 8.901 - 9.012; 
Os espaços em branco em torno dos sinais de =, +, *, ^, / e - são opcionais. Tais 
espaços são quase sempre usados para melhorar a legibilidade. Podem ser colocados 
diversos comandos na mesma linha desde que estejam separados por vírgulas ou ponto-e-
vírgulas. São exemplos: 
plot (x,y, ‘o’), text(1 ,20,’Sistema 1’), text(1 ,1 5,’Sistema 2’) 
plot(x,y, ‘o’); text(1 ,20, ‘Sistema 1’); text(1 ,1 5, ‘Sistema 2’) 
./ É o operador divisão elemento por elemento. B./A denota a divisão elemento por 
elemento, sendo requerido que A e B tenham a mesma dimensão. 
.* É o operador multiplicação elemento por elemento. X.*Y denota a multiplicação 
elemento por elemento, sendo requerido que A e B tenham a mesma dimensão. 
.^ É a potenciação elemento por elemento. Z = X.^Y denota a potenciação elemento 
por elemento, sendo requerido que X e Y tenham a mesma dimensão. 
O comando cd muda o diretório de trabalho corrente. Por exemplo, cd aula 
transforma o diretório corrente no especificado (aula). Por exemplo, cd .. acarreta a 
movimentação para o diretório acima do diretório corrente. O comando cd ele próprio 
informa o diretório corrente. 
O comando dir lista os diretórios. Por exemplo, dir aula lista os arquivos do diretório 
especificado (aula). 
O comando clc limpa a janela de comando e coloca o cursor no topo da tela. 
O comando help apresenta a documentação em linha do programa OCTAVE. O uso 
do comando help lista todos os tópicos primários, sendo que cada tópico primário 
corresponde a um nome de diretório do OCTAVE. O comando help tópico fornece auxílio no 
tópico especificado. O tópico pode ser um nome de comando ou um nome de diretório. Caso 
seja um comando, o help apresenta informações sobre o comando. Caso seja um diretório, o 
help apresenta a tabela de conteúdos do diretório. Não é necessário informar o caminho 
completo do diretório; sendo que o último componente, ou os últimos componentes são 
suficientes. Por exemplo, help general e help octave/general ambos listam a tabela de 
conteúdos do diretório toolbox/octave/general. 
O comando pause interrompe a execução de um programa até que o usuário digite 
alguma tecla, quando então o programa é executado a partir da primeira linha abaixo do 
comando pause. 
Par interromper definitivamente a execução de um programa devemos digitar ctrl c. 
13 
 
O OCTAVE possui um conjunto extenso de rotinas para obtenção de saídas gráficas. 
Se x e y forem vetores de mesmo comprimento, o comando plot(x,y) traça os valores 
de y contra os valores de x. Para traçar múltiplas curvas num único gráfico, utiliza-se o 
comando plot com argumentos múltiplos: 
plot(X1 ,Y1 ,X2,Y2,. . .,Xn,Yn) 
As variáveis X1 ,Y1 ,X2,Y2, e assim por diante, são pares de vetores. Cada par x-y é colocado 
no gráfico, resultando múltiplas curvas. Os argumentos múltiplos apresentam a vantagem de 
permitir que pares de vetores de diferentes comprimentos possam ser visualizados no 
mesmo gráfico. Cada par utiliza um tipo diferente de linha. O traçado de mais de uma curva 
num único gráfico também pode ser realizado utilizando-se o comando hold. O comando 
hold congela o gráfico atual. Portanto, as curvas subseqüentes serão sobrepostas à curva 
original. Digitando-se novamente o comando hold libera-se o gráfico presente. 
Uma vez que se tenha uma visualização gráfica na tela do monitor, é possível 
desenhar linhas de grade, colocar titulo no gráfico e colocar rótulos nos eixos x e y. Os 
comandos OCTAVE para grade, título, legenda do eixo dos x e legenda do eixo dos y são: 
grid 
title (‘título’) 
xlabel (‘legenda do eixo dos x’) 
ylabel (‘legenda do eixo dos y’) 
Para se escrever um texto começando no ponto (X, Y) da tela gráfica, se utiliza o 
comando: 
text(X, Y, ‘texto’) 
Por exemplo, o comando: text(3,0.45,’sen t’) escreverá horizontalmente a expressão sen t 
começando no ponto de coordenadas (3,0,45). Além disso, os comandos: 
plot(xl ,yl ,x2,y2), text(xl ,yl ,‘1‘), text(x2,y2,’2’) marcam as duas curvas de modo que elas 
possam ser facilmente distinguidas. 
plot(X,Y,’x’) traça um gráfico de pontos utilizando símbolos X, enquanto plot(X1 ,Y1 
,‘:‘,X2,Y2,’ +‘) utiliza uma linha de pontos (:) para a primeira curva e uma linha de símbolos 
(+) para a segunda curva. São os seguintes os tipos de linhas e de símbolos para as linhas de 
pontos disponíveis: 
Tipos de linha Tipos de ponto 
Cheia - Ponto . 
Tracejada -- Sinal mais + 
De pontos : Asterisco * 
Traço-ponto -. Circulo O 
 
plot(X,Y,’r’) 
plot(X,Y,’+g’) 
indicam o uso de uma linha vermelha no primeiro gráfico e de uma curva de marcas + em 
verde. As cores disponíveis são: 
 vermelho r 
 verde g 
 azul b 
 branco w 
 invisível i 
15 
 
 amarelo y 
 roxo m 
 azul-claro c 
No OCTAVE os gráficos são postos em escala automaticamente. Este gráfico 
permanece como gráfico corrente até que um outro seja traçado, situação em que o gráfico 
antigo é apagado e os eixos são automaticamente postos em nova escala. Em determinadas 
situações, pode ser desejável ignorar a característica de colocação automática dos eixos em 
escala, disponível no comando plot, e selecionar os pontos-limites manualmente. Para traçar 
uma curva numa região especificada por: 
v = [x-mín x-máx y-mín y-máx] 
deve ser digitado o comando axis(v). O comando axis(v), onde v é um vetor de quatro 
elementos, ajusta os eixos dentro dos limites especificados. Nos gráficos logarítmicos, os 
elementos de v são os 1og10 dos valores mínimo e máximo. A execução de axis(v) congela o 
corrente ajuste de escala dos eixos para os gráficos subseqüentes. Digitando-se axis 
novamente, retoma-se a colocação em escala automática. O comando axis(‘square’) faz com 
que o gráfico resultante fique numa tela em forma de quadrado. Com uma relação 
abscissa/ordenada de aspecto quadrado, uma retacom tangente igual a 1 é mostrada com 
inclinação verdadeira de 45°, e não distorcida pela forma irregular da tela do monitor. O 
comando axis(’normal’) faz a relação abscissa/ordenada original retornar. 
O comando subplot cria eixos em quadrantes específicos da janela de gráficos. 
subplot(m,n,p), ou subplot(mnp), quebra a janela gráfica, gerando matrizes m por n de eixos 
menores. Os eixos são contados (p) a partir da linha de topo da janela gráfica, e assim 
sucessivamente (segunda linha , etc). Por exemplo, 
subplot(2,1,1), plot(conv1) 
subplot(2,1,2), plot(conv2) 
Há 2 linhas e uma única coluna dividindo a janela gráfica, isto é, o gráfico da variável 
conv1 está na metade superior da janela gráfica e o gráfico da variável conv2 na metade 
inferior. 
Gráficos logarítmicos e polares são criados substituindo-se a expressão plot por 
loglog, semilogx, semilogy ou polar. Tais comandos são utilizados da mesma forma, 
modificando apenas como os eixos são postos em escala e como os dados são apresentados. 
O comando surf constrói uma superfície gráfica tridimensional. Por exemplo, 
surf(X,Y,Z,C) constrói o gráfico de uma superfície colorida cuja escala de cores é 
determinada pela faixa do argumento C. surf(X,Y,Z) emprega C = Z de modo que a cor é 
proporcional a altura da superfície. 
O comando whitebg transforma a cor de fundo da janela gráfica para branco. Os 
gráficos subseqüentes apresentarão também fundo branco. 
O comando bar constrói um gráfico de barras para os elementos do vetor Y: (bar(Y)). 
bar(X,Y) constrói um gráfico de barras dos elementos do vetor Y em localizações 
especificadas pelo vetor X. Os valores de X devem estar em ordem crescente. Caso os valores 
de X não estejam igualmente espaçados, o intervalo entre cada dado não é simétrico. As 
barras são colocadas no ponto médio entre dois valores adjacentes de X. 
O comando meshgrid gera arranjos X e Y para gráficos 3D. [X,Y] = meshgrid(x,y) 
transforma o domínio especificado pelos vetores x e y em arranjos X e Y que podem ser 
usados na avaliação de funções de duas variáveis e na construção de superfícies 3D. As 
linhas do arranjo X são cópias do vetor x e as colunas do arranjo Y são cópias do vetor y. 
17 
 
O comando plot3 traça linhas e pontos no espaço tridimensional. plot3() é o análogo 
tridimensional do plot(). plot3(x,y,z), onde x, y e z são três vetores de mesmo comprimento, 
traça uma linha no espaço 3-D através de pontos cujas coordenadas são os elementos de x, y 
e z. plot3(X,Y,Z), onde X, Y e Z são três matrizes de mesma dimensão, traça as diversas linhas 
das colunas de X, Y e Z. Vários tipos de linhas, de símbolos e de cores podem ser obtidos 
empregando plot3(X,Y,Z,s) onde s é análogo aos caracteres (cadeia alfanumérica) 
apresentados no comando plot. 
O comando mesh traça uma superfície reticulada tridimensional. mesh(X,Y,Z,C) traça 
uma rede 3-D (X,Y,Z) parametrizada por coloração (C). mesh(X,Y,Z) emprega C = Z, de forma 
que a coloração seja proporcional à altura da rede tridimensional. 
O comando colormap define as estruturas de cores, por exemplo: 
colormap(gray) - mapa de cores em escala linear de cinza. 
colormap(hot) - mapa de cores preto-vermelho-amarelo-branco. 
colormap(copper) - mapa de cores com tons lineares de cobre. 
colormap(prism) - mapa de cores do prisma. 
O OCTAVE possui um conjunto extenso de funções predefinidas e rotinas para a 
obtenção de saídas gráficas que podem ser chamadas pelo usuário para resolver diferentes 
tipos de problemas. 
2.2. Exemplos de programas no OCTAVE 
2.2.1. Programa comando2.m 
% Modelo para secagem de soja 
% Ys=umidade do sólido em um tempo t qualquer 
% Yse=umidade do sólido em tempo t -> infinito 
% Ys0=umidade do sólido em tempo t -> 0 
% Psat=pressão de saturação 
% Pv=pressão parcial de vapor no ar 
% t= tempo de secagem 
% m,n,q=parâmetros 
% Sistema: M=(Ys-Yse)/(Ys0-Yse)=exp[-m*((Psat-Pv)^n)*t^q] 
% Dados experimentais: 
% (Psat-Pv), (kgf/dm^2) 15 40 67 88 
% t, (h) 0,7 1,5 2,5 2,9 
% Sendo: M=(Ys-Yse)/(Ys0-Yse) 
m=0.1314; 
n=0.4601; 
q=0.0638; 
Psat_Pv=[15 40 67 88]; 
t=[0.7 1.5 2.5 2.9]; 
M=exp(-m*((Psat_Pv).^n).*t.^q) 
2.2.2. Programa comando3.m 
X1=[0 1 2 3 4 5 6 7 8 9 10] 
Y1=rand(1,11) 
X2=[0 1 2 3 4 5] 
Y2=rand(1,6) 
plot (X1,Y1, 'or'), text(4 ,0.8,'Sistema 1'), text(4 ,0.4,'Sistema 2') 
pause 
plot(X2,Y2, 'ob'); text(1 ,0.8, 'Sistema 1'); text(4 ,0.9, 'Sistema 2') 
pause 
subplot (2,2,1) 
plot(X1,'r'); 
subplot(2,2,2) 
plot(X2,'b'); 
subplot(2,2,3) 
plot(Y1,'g'); 
subplot(2,2,4) 
plot(Y2,'m'); 
pause 
clf 
plot(X1,Y1,X2,Y2), text(X1 ,Y1 ,'1'), text(X2,Y2,'2') 
pause 
grid 
title ('título') 
xlabel ('legenda do eixo dos x') 
ylabel ('legenda do eixo dos y') 
pause 
axis([0 5 0 1]); 
19 
 
2.2.3. Programa comando4.m 
whitebg 
X3=ones(1,3) 
Y3=[1 10 100] 
loglog(X3,Y3,'g') 
pause 
xx = [-5:0.25:5]; 
yy = [-5:0.25:5]; 
[x,y]=meshgrid(xx,yy) ; 
z = exp(-0.3*((x-2).^2.+(y-2).^2))+2*exp(-0.3*((x+1).^2.+(y-2).^2))'; 
surf(x,y,z) 
hold on 
xlabel('legenda do eixo dos x'); 
ylabel('legenda do eixo dos y'); 
zlabel('legenda do eixo dos z'); 
title('título'); 
pause 
clf 
X1=[0 1 2 3 4 5 6 7 8 9 10]; 
bar(X1) 
title ('título') 
xlabel ('legenda do eixo dos x') 
ylabel ('legenda do eixo dos y') 
pause 
figure 
plot(X1) 
title ('título') 
xlabel ('legenda do eixo dos x') 
ylabel ('legenda do eixo dos y') 
pause 
clf 
Vamos copiar a janela gráfica do OCTAVE para um documento do WORD. Vamos 
clicar em EDIT, depois COPY TO METAFILE, e em seguida, colamos no documento WORD, 
usando ctrl v. 
2.2.4. Programa comando5.m 
% Modelo para secagem de soja 
% Ys=umidade do sólido em um tempo t qualquer 
% Yse=umidade do sólido em tempo t -> infinito 
% Ys0=umidade do sólido em tempo t -> 0 
% Psat=pressão de saturação 
% Pv=pressão parcial de vapor no ar 
% t= tempo de secagem 
% m,n,q=parâmetros 
% Sistema: M=(Ys-Yse)/(Ys0-Yse)=exp[-m*((Psat-Pv)^n)*t^q] 
% Novos dados experimentais: 
% (Psat-Pv), (kgf/dm^2) 15 40 67 88 
% t, (h) 0,7 1,5 2,5 2,9 
% M .6398 0.4788 0.3811 0.3316 
% Sendo: M=(Ys-Yse)/(Ys0-Yse) 
Psat_Pv=[15 40 67 88]; 
t=[0.7 1.5 2.5 2.9]; 
M=[0.6398 0.4788 0.3811 0.3316]; 
plot3(t,Psat_Pv,M,'-.r') 
xlabel('Tempo [h]'); 
ylabel('Psat-Pv [kgf/dm^2]'); 
zlabel('(Ys-Yse)/(Ys0-Yse)'); 
title( [ 'Secagem de soja ' ] ); 
grid 
pause 
Psat_Pv=[15+15*rand(100,100) 40+40*rand(100,100) 67+67*rand(100,100) 
88+88*rand(100,100)]; 
t=[0.7+0.7*rand(100,100) 1.5+1.5*rand(100,100) 2.5+2.5*rand(100,100) 
2.9+2.9*rand(100,100)]; 
m=0.1314; 
n=0.4601; 
q=0.0638; 
M=exp(-m*((Psat_Pv).^n).*t.^q); 
colormap(prism); 
mesh(t,M,Psat_Pv) 
grid 
xlabel('Tempo [h]'); 
ylabel('(Ys-Yse)/(Ys0-Yse)'); 
zlabel('Psat-Pv [kgf/dm^2]'); 
title( [ 'Secagem de soja ' ] ); 
2.3. Exercício proposto 1 
 
 
 
 
 
21 
 
Olhar o help dos comandos e funções da Tabela 1: 
Na Tabela 1 são listados alguns destes comandos e funções. 
Tabela 1 
Comandos e Funções Matriciais Comumente 
Usados na Solução de Problemas 
Explicações sobre o que o Comando Faz, o que a 
Função Matricial Significa, ou o que a Declaração 
Significa 
ans 
atan 
axis 
clear 
clg 
conj 
cos 
cosh 
det 
diag 
 
exit 
exp 
eye 
format Iong 
format short 
format short e 
grid 
hold 
i 
imag 
inf 
inv 
Ienght 
linspace 
Iog 
Ioglog 
Iogspace 
log10 
max 
min 
NaN 
ones 
pi 
plot 
prod 
quit 
rand 
quit 
rand 
rank 
real 
rem 
sin 
 
sinh 
size 
sqrt 
sum 
tan 
tanh 
text 
title 
trace 
who 
xlabel 
ylabel 
zlabel 
zeros 
 
23 
 
2.4. Exercício proposto: Programa rpeq2.m 
Três correntes gasosas são alimentadas em um misturador. Sendo conhecidas as 
composições dessas três correntes e da corrente de saída, determinar as vazões molares de 
todas as correntes evolvidas no processo. 
 
Gás % molar corrente A 
% molar corrente 
B 
% molar corrente 
C 
% molar corrente 
D 
4CH 25% 35% 55% 40% 
62HC 35% 20% 40% 35%83HC 40% 45% 5% 25% 
 
misturador
A=100lbmol/h
B
C
D
 
Balanço global: DCBA  . 
Balanço por componente: 
D25,0C05,0B45,0A40,0
D35,0C40,0B20,0A35,0
D40,0C55,0B35,0A25,0



 
Análise: número de variáveis=3; número de equações=3 (o balanço global é a soma dos 
balanços por componente). 
3. Aula 3 – Introdução à programação em linguagem estruturada 
Queridos alunos e alunas, 
Tudo bem? 
3.1. Comandos 
No OCTAVE há seis operadores relacionais: <, <=, >, >=, ==, e ~=. A<B efetua 
comparações entre A e B, elemento por elemento, e retorna uma matriz de elementos 
unitários, caso a relação seja verdadeira, e elementos nulos, caso a relação seja falsa. A e B 
devem apresentar a mesma dimensão. 
& indica o e lógico. A & B é uma matriz cujos elementos são unitários, quando ambas 
as matrizes A e B apresentarem elementos não nulos; ou é uma matriz com elemento zero, 
quando tanto A quanto B apresentarem um elemento nulo. A e B devem possuir a mesma 
dimensão. 
| indica o ou lógico. A | B é uma matriz cujos elementos são unitários, quando A ou B 
apresentar um elemento não nulo, e nulos, quando ambas as matrizes A e B apresentarem 
elementos nulos. A e B devem apresentar a mesma dimensão a menos que sejam escalares. 
~ indica o complemento lógico não. ~A é uma matriz cujos elementos são unitários, 
quando A possuir elementos nulos, ou uma matriz cujos elementos são nulos, quando A 
possuir elementos não nulos. 
xor indica ou exclusivo. xor(A,B) é unitário quando A ou B, porém não 
simultaneamente, apresentarem elementos não nulos. 
 
 
25 
 
O comando while repete sentenças um número indefinido de vezes. A forma geral do 
comando while é: 
while variável 
Sentenças 
end 
As sentenças são executadas enquanto a variável apresentar todos os elementos não 
nulos. A variável é o resultado de: expressão ==, <, >, <=, >=, ~= expressão. Por exemplo: 
A=5; B= 1; 
while A > B, 
A=A-1; 
end 
O comando if executa sentenças condicionalmente. A forma genérica é: 
if variável 
Sentenças 
end 
As sentenças são executadas caso a parte real das variáveis apresente elementos não 
nulos. A variável é usualmente o resultado de: expressão ==, <, >, <=, >=, ~= expressão. Por 
exemplo: 
I=4; J=3; 
if I < J 
A(I,J) = 2; 
elseif abs(I-J) = = 1 
A(I,J) = -1; 
else 
A(I,J) = 0; 
End 
 
3.2. Exemplo de programas 
3.2.1. Programa operrel1.m 
m=10; 
w=20; 
m>w %errado (0) 
m<w %certo (1) 
a=11; 
while a > m & a < w 
a=a+1 
end 
m|w %certo (1) 
a=30 
while a > m | a > w 
a=a-1 
end 
~m %elemento diferente de zero (0) 
~ w %elemento diferente de zero (0) 
a=1; 
while a ~= m 
a=a+1 
end 
xor(m,w) %ambos são diferentes de zero (0) 
a=1; 
while xor(a>m,a<w) 
%a > m xor a < w 
a=a+1 
end 
3.2.2. Programa operrel2.m 
A=15; 
B=10; 
while A>B 
A=A-1 
end 
27 
 
pause 
while A>5 
A=A-1 
end 
pause 
while B>A 
B=B-1 
end 
clear all 
pause 
I=4 
J=3 
if I < J 
A(I,J)=2 
elseif abs(I-J) == 1 
A(I,J)=-1 
else 
A(I,J)=0 
end 
clear A 
pause 
I=3 
J=4 
if I < J 
A(I,J)=2 
elseif abs(I-J) == 1 
A(I,J)=-1 
else 
A(I,J)=0 
end 
3.3. Exercício proposto (aula): Programa operrel3.m 
Para as matrizes x e y escrever na tela do computador: 
 as matrizes x e y têm diagonais idênticas (x11=y11 e x22=y22); 
 as matrizes x e y têm diagonais distintas, porém, com um elemento idêntico; 
 as matrizes x e y têm diagonais distintas 
a) 











41
41
x , 










42
31
y 
b) 











41
41
x , 










42
31
y 
c) 











41
41
x , 










42
31
y 
 
29 
 
4. Aula 4 – Desenvolvimento de programas aplicados 
Queridos alunos e alunas, 
Tudo bem? 
Revisão sobre matrizes 
Inversa: 
A inversa de uma matriz A , definida apenas para matrizes quadradas, recebe a 
denominação 1A . 
  A
A
1
A 1 , sendo A a transposta da matriz cujos elementos ija de A são 
substituídos pelos respectivos cofatores ijA . 
A inversa de A somente existe se 0A  . 
A matriz quadrada A é singular se 0A  e não singular se 0A  . 
Exemplo: 









2221
1211
aa
aa
A ; 21122211 aaaaA  ; 
Vamos relembrar a definição de cofator, ijA : 
  ijjiij M1A  , sendo ijM o determinante da matriz resultante da eliminação da linha i 
e da coluna j da matriz A. 
  22221111 aa1A   
  21212112 aa1A   
  12121221 aa1A   
  11112222 aa1A   











1121
1222
aa
aa
A 











1121
1222
1
aa
aa
A
1
A 
Vamos analisar um exemplo numérico: 









32
01
A ; 










12
03
A ; 










3
1
3
2
01
A 1 
Posto: 
O posto ou rank de uma matriz é definido tanto para matrizes quadradas quanto para 
matrizes não quadradas. A maior sub-matriz quadrada com determinante não nulo fornece o 
posto da matriz. O posto é o número de linhas que é igual ao número de colunas desta 
submatriz. O posto de uma matriz tem muita utilidade na análise de sistemas de equações 
algébricas lineares: 
bAx  
Vamos definir a matriz aumentada do sistema como: 
 bAA,  
Ánálise: 
1) se r(A) diferente r(A'): o sistema linear não tem solução. 
2) se r(A) igual r(A'), então 
2a) se r(A)=r(A')<número de equações: o sistema linear tem infinitas 
soluções. 
2b) se r(A)=r(A')=número de equações: o sistema linear tem uma única 
solução. 
31 
 
4.1. Comandos do OCTAVE 
No OCTAVE, as funções: 
ones(n) 
ones(m,n) 
ones(A) 
geram matrizes especiais. Isto é, ones(n) produz uma matriz n X n de elementos iguais a um. 
ones(m,n) produz uma matriz m X n de elementos iguais a um. De modo semelhante, a 
função zeros(n) produz uma matriz n X n de zeros, enquanto zeros(m,n) produz uma matriz 
m X n de zeros. A função zeros(A) produz uma matriz de zeros com as mesmas dimensões de 
A. 
Se x é um vetor, a sentença diag(x) produz uma matriz diagonal com os elementos de 
x na diagonal principal. Por exemplo, para o vetor x=[ones(1,n)], diag([ones(1 ,n)]) fornece 
uma matriz identidade n X n. 
Se A for uma matriz quadrada, então diag(A) é um vetor formado pelos elementos da 
diagonal principal de A e diag(diag(A)) é uma matriz diagonal com elementos de diag(A) 
aparecendo na diagonal. 
O comando SUM indica a soma de elementos. Para vetores, SUM(X) é a soma de 
elementos de X. Para matrizes, SUM(X) é um vetor linha cujos elementos são a soma de cada 
coluna de X. O comando SUM(DIAG(X)) é o traço de X. 
O comando TRACE(A) é a soma dos elementos diagonais de A ou a soma dos 
autovalores de A. 
INV(A) é a inversa da matriz quadrada A. Uma mensagem de erro é impressa na tela 
se A for singular. 
RANK(A) é o número de linhas ou colunas linearmente independentes. 
DET(A) é o determinante da matriz quadrada A. 
Em programas OCTAVE é necessário freqüentemente entrar com a matriz identidade 
I. Uma sentença eye(n) fornece uma matriz identidade n X n. 
O formato STRINGS atribui a uma variável caracteres escritos dentro das aspas, por 
exemplo: X = 'Any Characters'. 
O comando INPUT fornece ao usuário o prompt e aguarda a entrada da variável 
digitada através do uso do teclado. 
No OCTAVE o comando global X Y Z define X, Y, Z em escala global. Tradicionalmente, 
cada função do OCTAVE, definida por um arquivo com extensão m apresenta suas próprias 
variáveis locais que são separadas das variáveis de outras funções ou do espaço de trabalho 
ou demais sub-rotinas. Entretanto, caso várias funções, e possivelmente o espaço de 
trabalho, declarem uma variável particular como global, então, todas compartilham uma 
cópia da variável. Qualquer atribuição a esta variável em qualquer função estará disponível 
para todas as demais funções que declararam a variável como global. 
O comando isglobal é verdadeiro para variáveis globais. ISGLOBAL(A) é unitário caso 
A seja uma variável global ou zero no casocontrário. 
Caso F seja uma linha contendo o nome de uma função, usualmente um arquivo.m, 
então FEVAL(F,x1,...,xn) calcula a função com os argumentos fornecidos entre parênteses. 
Por exemplo, F = 'foo', FEVAL(F,9.64) é o mesmo que foo(9.64). 
O comando NUM2STR efetua a transformação de um número em série de caracteres 
alfanuméricos. T = NUM2STR(X) converte o número escalar X em uma representação 
alfanumérica T com 4 dígitos e um expoente caso seja necessário. Este comando é útil na 
rotulagem de gráficos com os comandos TITLE, XLABEL, YLABEL, TEXT. 
33 
 
O comando FOR repete sentenças um número específico de vezes. A forma geral do 
FOR É: FOR variável = expressão, sentença, ..., sentença END. A expressão é freqüentemente 
da forma X:Y. Por exemplo, admitindo que N seja conhecido previamente: 
FOR I = 1:N, 
FOR J = 1:N, 
A(I,J) = 1/(I+J-1); 
END 
END 
FOR S = 1.0: -0.1: 0.0, END indica passos de S com incremento de -0.1. 
O comando SAVE salva as variáveis da área de trabalho no disco. O comando SAVE 
fname salva todas as variáveis da área de trabalho no arquivo binário "MAT-file" nomeado 
fname.mat. Os dados podem ser carregados com o comando LOAD. Caso haja omissão do 
nome do arquivo o comando SAVE utiliza o nome de arquivo padrão denominado 
"octave.mat". Por exemplo, 
SAVE fname X (salva apenas X). 
SAVE fname X Y Z (salva X, Y e Z). 
SAVE fname X Y Z -ascii (usa forma ASCII com 8-dígitos em vez de forma binária). 
SAVE fname X Y Z -ascii -double (usa forma ASCII com 16-dígitos). 
O comando LOAD restaura variáveis gravadas no disco. LOAD fname restaura as 
variáveis do arquivo MAT 'fname.mat'. O uso do comando LOAD carrega variável do arquivo 
'octave.mat'. LOAD xxx.yyy lê o arquivo ASCII xxx.yyy, que contém um arranjo retangular de 
dados numéricos, organizados em m linhas e n valores em cada linha, sendo o resultado uma 
matriz m por n denominada xxx. 
4.2. Exemplo de programas 
4.2.1. Programa comand61.m 
n=3; 
y=[zeros(1,n)] 
x=[ones(1,n)] 
diag([ones(1 ,n)]) 
pause 
a=[10 2 3;4 59 6;57 86 29] 
diag(a) 
diag(diag(a)) 
trace(a) 
sum(diag(a)) 
rank(a) 
det(a) 
inv(a) 
4.2.2. Programa Comand62.m 
Avaliar a existência e unicidade dos sistemas: 
a) 




7y2x2
3yx ; b) 




6y2x2
3yx ; c) 




4yx2
3yx ; d) 




0yx2
0yx ; e) 




0y2x2
0yx 
Programa: 
'informar o número de equações' 
n=input('O número de equações é ====>'); 
'informar a matriz de coeficientes' 
mat = input ('A matriz de coeficientes é ====> '); 
r1=rank(mat); 
'informar o vetor de constantes' 
x = input ('O vetor de constantes é ====> '); 
mata=[mat x] 
r2=rank(mata); 
if r1 ~= r2 
'O sistema linear não tem solução' 
'O rank da matriz de coeficiente é diferente' 
 'do rank da matriz aumentada' 
end 
% para o caso de r1 == r2 
 if r1 == r2 
 if r2 < n 
 'Há infinitas soluções' 
'O rank da matriz de coeficiente ou o rank da matriz aumentada' 
35 
 
'é menor que o número de equações' 
end 
 if r2 == n 
 'Há uma única solução' 
'O rank da matriz de coeficiente ou o rank da matriz aumentada' 
'é igual ao número de equações' 
 x=mat\x 
end 
end 
4.3. Exercício proposto: comand63.m 
A matriz atômica (A) indica no elemento ija o número de átomos do elemento 
químico i presente na molécula de um composto j. Dado um conjunto de NC compostos, o 
número de reações químicas independentes (NR) que podem ser geradas é igual a: 
)A(postoNCNR  
a) Determinar o número de reações independentes para o seguinte conjunto de 
compostos: 
OH;ClHC;ClHC;O;HCl;Cl;HC 2322422242 
Os coeficientes estequiométricos das reações químicas independentes podem ser 
calculados pela seguinte equação matricial: 
 





































47
37
27
17
46
36
26
16
45
35
25
15
43
33
23
13
42
32
22
12
41
31
21
11
44434241
34333231
24232221
14131211
a
a
a
a
a
a
a
a
a
a
a
a
v
v
v
v
v
v
v
v
v
v
v
v
aaaa
aaaa
aaaa
aaaa
 
 
onde ija denotam os elementos da matriz atômica A e ijv representam os 
coeficientes estequiométricos do composto i na reação k. Os coeficientes estequiométricos 
dos demais compostos são dados por: 
 





















100
010
001
vvv
vvv
vvv
737271
636261
535251
 
 
b) Apresentar um conjunto de reações químicas independentes com os seus 
respectivos coeficientes estequiométricos. 
37 
 
5. Aula 5 - Métodos computacionais aplicados à resolução de sistemas lineares de 
equações algébricas 
 Queridos alunos e alunas, 
Tudo bem? 
5.1. Comandos 
O comando GTEXT posiciona um texto em um gráfico 2-D usando o mouse do 
computador. Por exemplo, GTEXT ('caractere') requer o acionamento de uma letra do 
teclado ou do botão do mouse. 
5.2. Método Gauss-Seidel 
Resolução de sistemas lineares de equações algébricas - Método Gauss-Seidel. O 
método Gauss-Seidel é empregado em sistemas predominantemente diagonais. Cada 
equação é usada para calcular a variável da diagonal. É um método iterativo que parte de 
uma estimativa inicial. Sempre devemos usar os valores mais recentes de cada variável. 
 
Exemplo1: Reatores CSTR em série 
 
 
V1 V2 V3 V4
CA1 CA2 CA3 CA4
CA0=1 mol/l
1000 l/h
1000 l/h
1100 l/h
100 l/h
1100 l/h
100 l/h
1100 l/h 1000 l/h
 
 
100 l/h 
Hipóteses: 
Estado estacionário. 
As densidades de todas as correntes são iguais, então ao invés de balanço de massa 
fazer balanço de vazão. 
Todas as correntes são líquidas. 
Reação: A=>B (reação de primeira ordem) 
A dependência entre a taxa e a concentração é linear 
Taxa de consumo de A, sendo i o reator i: Ri=Ki*CAi*Vi 
Sistema linear de equações resultante é: 














01250C-1100C
0100C1240C-1100C
0100C1400C-1000C
1000C1001000C
A4A3
A4A3A2
A32A1A
1A1A
 
Dados cinéticos: 
Reator Vi(L) Ki(h-1) 
1 1000 0,1 
2 1500 0,2 
3 100 0,4 
4 500 0,3 
 
O sistema é predominantemente diagonal porque os reciclos são pequenos quando 
comparados com o fluxo principal: 
39 
 
1100> 0 
-1400>1000+100 
-1240>1100+100 
-1250>1100 
5.3. Exemplo de programa 
5.3.1. Programa gseidelq.m 
 
% Exemplo: i Reatores CSTR em série 
% Hipótes: 
% Estado estacionário. 
% As densidades de todas as correntes são iguais, então 
% ao invés de balanço de massa fazer balanço de vazão. 
% Todas as correntes são líquidas. 
% Reação: A=>B (reação de primeira ordem) 
% A dependência entre a taxa e a concentração é linear 
% Taxa de consumo de A, sendo i o reator i: 
% Ri=Ki*CAi*Vi 
% Sistema linear de equações resultante é: 
% 1100*CA1=1000 
% 1000*CA1-1400*CA2+100*CA3=0 
% 1100*CA2-1240*CA3+100*CA4=0 
% 1100*CA3-1250*CA4=0 
% O sistema é predominantemente diagonal porque os reciclos são pequenos quando 
comparados com o fluxo principal: 
% 1100> 0 
% -1400>1000+100 
% -1240>1100+100 
% -1250>1100 
%****************************************************************** 
% Como CA1=1000/1100, então, CA2, CA3, CA4 devem ser menor que CA1 
% CA2, CA3, CA4 devem estar no intervalo (0,CA1) 
% Estimativa inicial primitiva: CA2=CA3=CA4=0,6 
% Chute inicial 
CA_2=0.6; 
CA_3=0.6; 
CA_4=0.6; 
CA2old=CA_2; 
CA3old=CA_3; 
CA4old=CA_4; 
CA2=0.6; 
CA3=0.6; 
CA4=0.6; 
CA1=1000/1100; 
CA2=(1000*CA1+100*CA3)/1400; 
CA3=(1100*CA2+100*CA4)/1240; 
CA4=1100*CA3/1250; 
% Dado inicial 
CA_1=1000/1100; 
k=1; 
while abs((CA2-CA2old)/CA2) > 1e-7 | abs((CA3-CA3old)/CA3) > 1e-7 | abs((CA4-
CA4old)/CA4) > 1e-7 
k 
% Cada equação é usada para calcular a variável da diagonal: 
CA1=1000/1100; 
CA2=(1000*CA1+100*CA3)/1400; 
CA3=(1100*CA2+100*CA4)/1240; 
CA4=1100*CA3/1250; 
CA_1=[CA_1,CA1]; 
CA_2=[CA_2,CA2]; 
CA_3=[CA_3,CA3]; 
CA_4=[CA_4,CA4]; 
l=length(CA_1)-1; 
CA2old=CA_2(l); 
CA3old=CA_3(l); 
CA4old=CA_4(l); 
k=k+1; 
end 
CA=[CA_1' CA_2' CA_3' CA_4']; 
% saídas gráficas 
subplot (2,2,1) 
plot(CA(:,1),'r');xlabel('Iteração'); 
ylabel('CA1'); 
subplot(2,2,2) 
plot(CA(:,2),'b'); 
xlabel('Iteração'); 
ylabel('CA2'); 
subplot(2,2,3) 
plot(CA(:,3),'g'); 
xlabel('Iteração'); 
ylabel('CA3'); 
subplot(2,2,4) 
plot(CA(:,4),'m'); 
41 
 
xlabel('Iteração'); 
ylabel('CA4'); 
gtext('CSTR - Método de Gauss-Seidel'); 
% salvar dados 
save gseidel.dat CA -ascii 
5.4. Exercício proposto: Programa gseidelt.m 
A determinação do perfil de temperatura de um sólido no estado estacionário pode 
ser expressa por um sistema de equações algébricas lineares. Para calcular o perfil de 
temperatura do corpo bidimensional devem ser resolvidas as seguintes equações algébricas 
lineares, sendo T1, T2, ... Tn as temperaturas pontuais a serem determinadas: 
-4*T1+T2+T4=-600 
T1-4*T2+T3+T5=-500 
2*T2-4.67*T3+T6=-567 
T1-4*T4+T5+T7=-100 
T2+T4-5*T5+T6+T8=0 
T3+2*T5-4.67*T6+T9=-67 
T4-4.67*T7+T8=-167 
2*T5+T7-4.67*T8+T9=-67 
T8-2.67*T9=-67 
 
6. Aula 6 – Métodos computacionais aplicados à resolução de sistemas lineares de 
equações algébricas 
Queridos alunos e alunas, 
Tudo bem? 
Vamos estudar o Método de Eliminação Gaussiana 
6.1. Método de Eliminação Gaussiana 
 A técnica de Eliminação Gaussiana consiste em realizar operações elementares 
partindo de bAx  : 





































n
2
1
n
2
1
nn2n1n
n22221
n11211
b
b
b
x
x
x
aaa
aaa
aaa





 até chegar à forma: 











































,
n
,
2
,
1
n
1n
2
1
nn
n,1n1n,1n
n222
n11211
bb
bb
bb
x
x
x
x
aa000
aaaa
aaaa0
aaaaaa




. Na forma triangular é fácil resolver o sistema 
por retrossubstituição: 
 na última linha, determinamos nx ; 
 conhecendo nx , na penúltima linha, calculamos 1nx  ; 
 conhecendo nx , 1nx  , calculamos 2nx  na linha anterior; 
 etc..., até determinar 1x na primeira linha. 
6.2. Comandos 
\ É o operador divisão esquerda. A\B equivale a fazer INV(A)*B. Caso A seja uma 
matriz N por N e B um vetor coluna com N componentes, então, X = A\B é a solução da 
equação A*X = B calculada por Eliminação Gaussiana. A\EYE(SIZE(A)) produz a inversa de A. 
43 
 
 
6.3. Exercício proposto: eg.m 
Cinco correntes gasosas são alimentadas em um misturador. Sendo conhecidas as 
composições dessas cinco correntes e da corrente de saída, determinar as vazões molares de 
todas as correntes evolvidas no processo. 
Resolver implementando \ e o método da Eliminação Gaussiana. 
 
misturador
A=100lbmol/h
B C
D
E F 
 
Gás % molar corrente A 
% molar 
corrente B 
% molar 
corrente C 
% molar 
corrente D 
% molar 
corrente 
E 
% molar 
corrente 
F 
4CH 39% 30% 40% 32% 28% 35% 
62HC 28% 25% 30% 23% 25% 20% 
83HC 23% 15% 10% 20% 20% 23% 
2H 9% 10% 11% 14% 16% 15% 
2N 1% 20% 9% 11% 11% 7% 
7. Aula 7 – Métodos computacionais aplicados à resolução de sistemas lineares de 
equações algébricas 
Queridos alunos e alunas, 
Tudo bem? 
Vamos estudar o Método Gauss-Jordan 
7.1. Método Gauss-Jordan 
 Para o sistema bAx  , após uma série de operações elementares: cIx  . 
   1selementareoperações AsoluçãoIIbA   
 Este método permite calcular a matriz inversa. Mas a matriz inversa não é sempre 
necessária. Se 1A  não é necessária, é melhor não usar o método de Gauss-Jordan porque é 
possível resolver o sistema linear por outros métodos (por exemplo, Eliminação Gaussiana) 
com menos operações algébricas 
7.2. Exercício proposto (aula): gj1.m 
 Uma bateria de três extratores opera em contra-corrente a fim de recuperar o produto 
P que está dissolvido no solvente A. Um solvente B é utilizado para a extração. Devido a um 
problema operacional no extrator 2, ele opera com capacidade reduzida, de modo que foi 
montado um bypass parcial que desvia 20% da vazão de saída do extrator 1. Admite-se que 
os solventes A e B sejam completamente imiscíveis e que o processo opere em estado 
estacionário. Calcular as concentrações de P em todas as etapas do processo, sabendo que 
as vazões dos dois solventes Vb,Va são iguais, que as concentrações de P na entrada da 
bateria são de: 
BsolventedekgPdekg01,04y  
45 
 
AsolventedekgPdekg2,00x  
 A equação de equilíbrio que relaciona as concentrações de saída de cada extrator é: 
x5,1y  
 
Extrator
1
Extrator
2
Extrator
3
Vb, y1
Va, x0
y2
x1 x1
20% Va, x1
x2
y3 y4, Vb
x3, Va
 
7.3. Exercício proposto (aula): gj2.m 
 A relação entre a pressão de vapor satP de um fluido e sua temperatura pode ser 
modelada pela equação: 


  TlnC
T
B
AexpP sat , na qual T é expressa em Rankine e satP 
em psia. Na tabela abaixo são fornecidos dados experimentais de T versus satP para a 
amônia. 
 
T(°R) 459,67 479,67 499,67 
satP (psia) 30,42 48,21 73,32 
 
 Calcular a pressão de vapor da amônia a 474,67 °R. 
 
8. Aula 8 – Ajuste e interpolação de curvas oriundas de resultados experimentais 
Queridos alunos e alunas, 
Tudo bem? 
Ajuste 
Qual o conjunto de parâmetros que fornece a melhor representação dos dados 
experimentais segundo um modelo previamente escolhido? 
Ajuste por regressão linear 
O objetivo é minimizar:  


NP
1j
2
jj yYf , sendo 
NP =número de pontos experimentais; 
jY =valor experimental da variável dependente no ponto j; 
jy =valor calculado da variável dependente no ponto j; 
Vamos propor um modelo linear nos parâmetros n10 b,,b,b  : 
     jnnj11j00j xgbxgbxgby   , sendo: 
jx =valor da variável independente no ponto experimental j 
Vamos substituir na função f o valor de jy : 
            
 







NP
1j
2n
0i
jiij
NP
1j
2
jnnj11j00j
NP
1j
2
jj xgbYxgbxgbxgbYyYf  
Vamos derivar f em relação aos parâmetros n10 b,,b,b  e igualar as derivadas a zero. 
Dessas operações resulta um sistema de equações algébricas lineares nos parâmetros 
n10 b,,b,b  com o formato:    YGbGG TT  , sendo 
47 
 
   
   
   


























NP
2
1
NPnNP1
2n21
1n11
Y
Y
Y
Y
xgxg1
xgxg1
xgxg1
G





 
Este sistema pode ser resolvido por eliminação Gaussiana. 
8.1. Comandos 
O comando POLYFIT(x,y,n) encontra os coeficientes do polinômio p(x) de grau n que 
se ajusta aos dados, p(x(i)) ~= y(i), por mínimos quadrados. [p] = POLYFIT(x,y,n) fornece os 
coeficientes polinomiais p. 
O comando LENGTH indica o número de componentes de um vetor. Por exemplo, 
LENGTH(X) retorna o comprimento do vetor X. Este comando é equivalente a MAX(SIZE(X)). 
8.2. Exemplo de programa 
8.2.1. Programa ajuste1.m 
Os seguintes dados experimentais estão disponíveis para a condutividade térmica da 
água à pressão de 50 bar: 
 
 
T (°C) 0 50 100 150 200 
k 
(mwatt/m°C) 573 647 684 690 668 
 
Ajustar um polinômio cúbico por regressão linear a partir de todas as informações 
disponíveis. 
% Aplicação da técnica de regressão linear ao ajuste de um 
%modelo linear nos parâmetros 
% Ajuste por regressão linear: 
% min f=somatório(Yj-yj)^2 
% Propondo-se um modelo linear: 
% yj=b0g0(xj)+b1g1(xj)+...+bngn(xj) 
% O conjunto de parâmetros é solução do sistema linear: 
% (G'G)b=G'Y - Eliminação Gaussiana 
% Sistema: Ajuste da condutividade térmica da água à pressão de 50 bar 
% k=b0+b1*T+b2*T^2+b3*T^3, sendo 
% k=condutividade térmica 
% T=temperatura 
% Dados experimentais: 
% T(°C) 0 50 100 150 200 
% k(mwatt/m°C) 573 647 684 690 668 
T=[0 50 100 150 200]'; 
k=[573 647 684 690 668]'; 
[b]=polyfit(T,k,3) 
% Alternativamente 
% Modelo linear:k=b0+b1*T+b2*T^2+b3*T^3 ou 
% Y=b0*g0+b1*g1+b2*g2+b3*g3 
% Sendo: Y=k; 
% g0=ones(1,size(k)); 
% g1=T; 
% g2=T.^2; 
% g3=T.^3; 
Y=k 
g0=ones(1,length(k))' 
g1=T 
g2=T.^2 
g3=T.^3 
G=[g0 g1 g2 g3] 
% (G'G)b=G'Y - Eliminação Gaussiana 
% Sendo: GTG=(G'G) e GTY=G'Y 
GTG=G'*GGTY=G'*Y 
% Conjunto de parâmetros que fornece a melhor representação 
% dos dados experimentais segundo um modelo previamente escolhido 
b=GTG\GTY 
flipud(b)' 
49 
 
8.3. Exercício proposto (aula): Programa ajuste2.m 
Sistema: Ajuste do custo de mão-de-obra na produção de trocadores de calor: 
 C=b0+b1*A+b2*N, sendo 
 C=custo da mão-de-obra 
 A=área do trocador 
 N=número de tubos 
 Dados experimentais: 
 C 310 300 275 250 220 200 190 150 140 100 
 A 120 130 108 110 84 90 80 55 64 50 
 N 550 600 520 420 400 300 230 120 190 100 
8.3.1. Programa ajuste3.m 
Ajuste de um modelo não linear nos parâmetros - Secagem de soja 
Modelo para secagem de soja 
Ys=umidade do sólido em um tempo t qualquer 
Yse=umidade do sólido em tempo t -> infinito 
Ys0=umidade do sólido em tempo t -> 0 
Psat=pressão de saturação 
Pv=pressão parcial de vapor no ar 
t= tempo de secagem 
m,n,q=parâmetros a ajustar 
Sistema: M=(Ys-Yse)/(Ys0-Yse)=exp[-m*((Psat-Pv)^n)*t^q] 
Dados experimentais: 
M, (adimensional) 0,674 0,561 0,405 0,327 
(Psat-Pv), (kgf/dm^2) 12 25 60 90 
t, (h) 0,5 1 2 3 
51 
 
9. Aula 9 – Ajuste e interpolação de curvas oriundas de resultados experimentais 
Queridos alunos e alunas, 
Tudo bem? 
Interpolação 
Na interpolação, a curva contínua obtida passa por todos os pontos experimentas. A 
interpolação é empregada quando não se tem um modelo para um dado fenômeno ou 
quando o modelo é complicado demais para as necessidades. Além disso, a interpolação 
polinomial tem importância na interpolação de dados tabelados. Um polinômio interpolador 
de grau n passa por n+1 pontos distintos. Este polinômio interpolador é único, embora haja 
maneiras diferentes de obtê-lo (resolução do sistema linear, forma de Lagrange, forma de 
Newton para pontos regularmente espaçados). Na prática raramente se usam polinômios 
interpoladores de grau muito elevado (>7 ou 8), pois o comportamento destes polinômios 
pode ser errático. 
Seja o polinômio interpolador de grau n: nn2210n xaxaxaaP   . Supõe-se 
que há n+1 pontos experimentais disponíveis:         ;xf,x;;xf,x;xf,x nn1100  e n+1 
coeficientes a determinar: n10 a;;a;a  . Vamos impor que o polinômio passe exatamente 
pelos pontos tabelados: 
 
 
 
 









n
n
nn
2
n2n10
1
n
1n
2
12110
0
n
0n
2
02010
xfxaxaxaa
xfxaxaxaa
xfxaxaxaa




 
 
n10 x;x;x  são conhecidos. Portanto, o problema é resolver um sistema de equações 
lineares em n10 a;;a;a  de forma matricial: 
 
 
 
 






































n
1
0
n
1
0
n
n
n
1
n
0
2
nn
2
11
2
00
xf
xf
xf
a
a
a
x
x
x
xx1
xx1
xx1





 
 
A matriz de coeficientes é denominada matriz de Vandermonde. Para n10 x;x;x  
todos mutuamente distintos, o posto da matriz é n+1, o que garante a unicidade da solução 
(o polinômio interpolador é único). 
Forma de Lagrange do Polinômio Interpolador 
             xLxfxLxfxLxfxp nn1100n  
Sendo:     
    
       nj1jj1jj1j0j
n1j1j10
j xxxxxxxxxx
xxxxxxxxxx
xL







 
Para   1xLxx jjj  ; 
Para 0Lx,,x,x,,x,xx jn1j1j10    . 
Compactando a notação:     

 


n
ji
0i ij
i
j xx
xx
xL e o polinômio é:      xLxfxp j
n
0j
jn 

 . 
9.1. Exemplo de programa 
9.1.1. Programa intp1.m 
Interpolação de dados em uma tabela de vapor. 
De acordo com os dados da tabela teremos um polinômio interpolador de grau 3 
passando por 4 pontos distintos:   3322103 tatataatP  . 
53 
 
 
373
T
t  
 





































462
1406
1939
2257
a
a
a
a
1228,59717,27239,11
6252,33599,25362,11
0392,26081,12681,11
1111
3
2
1
0
 
 
T(K) t (T adimensional) 
t2(T 
adimensional)2 
t3(T 
adimensional)3 
vapH 
373 1 1 1 2257 
473 1,2681 1,6081 2,0392 1939 
573 1,5362 2,3599 3,6252 1406 
643 1,7239 2,9717 5,1228 462 
 
% Interpolação de dados em uma tabela de vapor 
% Dados 
global a0 a1 a2 a3 
temp=[373 
 473 
 573 
 643]; 
DH_VAP=[2257 
 1939 
 1406 
 462]; 
% Polinômio interpolador de grau 3 a partir de 4 dados 
% P3_t=a0+a1*t+a2*t^2+a3*t^3 
A=[ones(4,1) temp/373 (temp/373).^2 (temp/373).^3] 
B=DH_VAP; 
a=A\B 
a0=a(1); 
a1=a(2); 
a2=a(3); 
a3=a(4); 
% Para T=512 K 
T=512; 
DH_512calc=a(1)+a(2)*T/373+a(3)*(T/373)^2+a(4)*(T/373)^3 
% Na tabela de vapor DH_512exp=1768 KJ/Kg 
DH_512exp=1768 
Desvio_relativoH=(DH_512calc-DH_512exp)/DH_512exp 
% Vamos extrapolar DH_VAP para o ponto crítico da água, onde DH_VAP=0, isto é, fase 
vapor e 
% líquida são identicas. Vale ressaltar que na extrapolação não há garantia do 
comportamento da função. 
DH_VAP_TC=0; 
TC_calc=fzero('intp11',650) 
% Dado experimental: TC_exp=647 K 
TC_exp=647; 
55 
 
Desvio_relativoTC=(TC_calc-TC_exp)/TC_exp 
% Erro na extrapolação maior que na interpolação 
plot(temp,DH_VAP,'ro'); 
hold 
plot(T,DH_512calc,'w+'); 
axis([300 900 0 2500]); 
text(550,DH_512calc,'Ponto interpolado'); 
hold 
hold 
plot(TC_calc,DH_VAP_TC,'w*'); 
v = [350 900 0 2500]; 
axis(v) 
text(TC_calc,150,'Ponto extrapolado'); 
title ('Interpolação/extrapolação de dados em tabela de vapor'); 
xlabel ('Temperatura (K)'); 
ylabel ('DH, KJ/Kg'); 
%**************************************************************************
**************************************** 
função intp11.m 
function [f]=intp11(x); 
 
global a0 a1 a2 a3 
f=a0+a1*x/373+a2*(x/373)^2+a3*(x/373)^3; 
9.1.2. Programa intp2.m 
% Interpolação de dados em uma tabela de vapor 
% Dados 
temp=[373 
 473 
 573 
 643]; 
DH_VAP=[2257 
 1939 
 1406 
 462]; 
% Polinômio interpolador de grau 3 a partir de 4 dados 
% P3_t=a0+a1*t+a2*t^2+a3*t^3 
% Para T=512 K 
T=512; 
DH_512calc=interp1(temp,DH_VAP,T,'spline') 
% Na tabela de vapor DH_512exp=1768 KJ/Kg 
DH_512exp=1768 
Desvio_relativoH=(DH_512calc-DH_512exp)/DH_512exp 
plot(temp,DH_VAP,'ro'); 
hold 
plot(T,DH_512calc,'w+'); 
axis([300 900 0 2500]); 
text(550,DH_512calc,'Ponto interpolado'); 
57 
 
title ('Interpolação/extrapolação de dados em tabela de vapor'); 
xlabel ('Temperatura (K)'); 
ylabel ('DH, KJ/Kg'); 
 
9.1.3. Exercício proposto: Programa lagrange.m (Dica: usar conv2) 
Interpolação de dados em uma tabela de vapor. 
De acordo com os dados da tabela teremos um polinômio interpolador de grau 3 
passando por 4 pontos distintos:   3322103 tatataatP  . 
373
T
t  
T(K) 
373
T
t  vapH 
373 1 2257 
473 1,2681 1939 
573 1,5362 1406 
643 1,7239 462 
 
         
               tLtHtLtHtLtHtLtH
tLtHtLtftp
33
vap
22
vap
11
vap
00
vap
j
3
0j
j
vap
j
3
0j
j3

 
 
        302010
321
0 tttttt
tttttt
tL


 ;         312101
320
1 tttttt
tttttt
tL


 ; 
        321202
310
2 tttttt
tttttt
tL


 ;         231303
210
3 tttttt
tttttt
tL


 
Programa: Determinar os coeficientes do polinômio e calcular  512H vap e comparar com o 
valor obtido da tabela de vapor:   kgkJ1768512H vap  . 
10. Aula 10 – Resolução de equação algébrica não linear 
Queridos alunos e alunas, 
Tudo bem? 
Método de Newton-Raphson 
O método se baseia em realizar linearizações sucessivas de  xf . Usando série de 
Taylor truncada no termo linear:      
0x
00 dx
df
xxxfxf  . Vamos impor que   0xf  , logo: 
      0xfxxxf 0,00  . Portanto, 
 
 0,
0
0
xf
xf
xx  . Generalizando,  
 i,
i
i1i
xf
xf
xx  . 
 
f(x)
xx0x1x2
tangente à curva 
em (x0, f(x0))
tangente à curva 
em (x1, f(x1))
tangente à curva 
em (x2, f(x2))
 
 
Críticado método de Newton-Raphson: 
 muito rápido desde que se tenha uma boa estimativa inicial 
 a convergência não é garantida 
 exige derivadas analíticas 
59 
 
10.1. Comandos 
O comando global X Y Z define X, Y, Z em escala global. Tradicionalmente, cada 
função, definida por um arquivo com extensão m apresenta suas próprias variáveis locais 
que são separadas das variáveis de outras funções ou do espaço de trabalho ou demais sub-
rotinas. Entretanto, caso várias funções, e possivelmente o espaço de trabalho, declarem 
uma variável particular como global, então, todas compartilham uma cópia da variável. 
Qualquer atribuição a esta variável em qualquer função estará disponível para todas as 
demais funções que declararam a variável como global. 
O comando isglobal é verdadeiro para variáveis globais. ISGLOBAL(A) é unitário caso 
A seja uma variável global ou zero no caso contrário. 
Caso F seja uma linha contendo o nome de uma função, usualmente um arquivo.m, 
então FEVAL(F,x1,...,xn) calcula a função com os argumentos fornecidos entre parênteses. 
Por exemplo, F = 'foo', FEVAL(F,9.64) é o mesmo que foo(9.64). 
O comando NUM2STR efetua a transformação de um número em série de caracteres 
alfanuméricos. T = NUM2STR(X) converte o número escalar X em uma representação 
alfanumérica T com 4 dígitos e um expoente caso seja necessário. Este comando é útil na 
rotulagem de gráficos com os comandos TITLE, XLABEL, YLABEL, TEXT. 
O comando FOR repete sentenças um número específico de vezes. A forma geral do 
FOR É: FOR variável = expressão, sentença, ..., sentença END. A expressão é freqüentemente 
da forma X:Y. Por exemplo, admitindo que N seja conhecido previamente: 
FOR I = 1:N, 
FOR J = 1:N, 
A(I,J) = 1/(I+J-1); 
END 
END 
FOR S = 1.0: -0.1: 0.0, END indica passos de S com incremento de 0.1. 
O comando SAVE salva as variáveis da área de trabalho no disco. O comando SAVE 
fname salva todas as variáveis da área de trabalho no arquivo binário "MAT-file" nomeado 
fname.mat. Os dados podem ser carregados com o comando LOAD. Caso haja omissão do 
nome do arquivo o comando SAVE utiliza o nome de arquivo padrão denominado "*.mat". 
Por exemplo, 
SAVE fname X (salva apenas X). 
SAVE fname X Y Z (salva X, Y e Z). 
SAVE fname X Y Z -ascii (usa forma ASCII com 8-digitos em vez de forma binária). 
SAVE fname X Y Z -ascii -double (usa forma ASCII com 16-digitos). 
O comando LOAD restaura variáveis gravadas no disco. LOAD fname restaura as 
variáveis do arquivo MAT 'fname.mat'. O uso do comando LOAD carrega variável do arquivo 
'*.mat'. LOAD xxx.yyy lê o arquivo ASCII xxx.yyy, que contém um arranjo retangular de dados 
numéricos, organizados em m linhas e n valores em cada linha, sendo o resultado uma 
matriz m por n denominada xxx. 
10.2. Métodos matemáticos para resolução de equações algébricas não lineares 
 
T, P
F, z
V, y
L, x 
61 
 
 
Flash isotérmico ( psia1600P  e F120T o ) 
Balanço por componente: iii VyLxFz  
Balanço global: VLF  
Restrições: 1z
i
i  , 1x
i
i  , 1y
i
i  
Relação de equilíbrio: iii xky  ,  P,Tkk ii  
Vamos dividir o balanço por componente por F : iii yF
V
x
F
L
z  
Sendo VLF  , 
F
V
1
F
L
 
Vamos definir  como a fração vaporizada da carga (alimentação), então: 
  iii yx1z  
Vamos usar a relação de equilíbrio, então: 
   iii xk1z   
i
i
i k1
z
x

 e 
i
ii
iii k1
zk
xky

 
Vamos usar 1x
i
i  e 1y
i
i  da seguinte forma: 
 
  01k1
z1k
k1
zzk
0xy
i i
ii
i i
iii
i
i
i
i  




 
Portanto a equação de Flash é:      01k1
z1k
f
i i
ii 


 
 
No do 
componente 
Componente iz ik 
1 metano 0,7 3,09 
2 propano 0,1 0,39 
3 n-pentano 0,1 0,093 
4 n-hexano 0,1 0,065 
 
10.3. Exemplo: Resolução de equações algébricas não lineares 
10.3.1. Programa new.m 
% Método de Newton-Raphson 
% Flash isotérmico; p=1600 psia; T=120°F 
% N Componentes: metano (zM=0,70); propano (zP=0,10); n-pentano (znP=0,10); n-
hexano (znH=0,10) 
% Constante de equilíbrio: metano (kM=3,09); propano (kP=0,39); n-pentano 
(knP=0,093); n-hexano (knH=0,065) 
% Equação de Flash: fb(i)=somatório[(ki-1)zi/(1+b(ki-1))]=0 
% Recursão do método: b(i+1)=b(i)-fb(i)/f'b(i) 
% f'b(i)=somatório[-(zi*(ki-1)^2)/(1+b*(ki-1))^2; 
% Estimativas iniciais: b=0,50 
global N 
N=4; 
z(1)=0.70; 
z(2)=0.10; 
z(3)=0.10; 
z(4)=0.10; 
k(1)=3.09; 
k(2)=0.39; 
k(3)=0.093; 
k(4)=0.065; 
b=0.50; 
b0=b; 
f=feval('new1',b,z,k); 
f0=f; 
df=feval('new2',b,z,k); 
b=b-(f/df); 
while abs(f) > 1e-15 
 b0=[b0,b]; 
 f=feval('new1',b,z,k); 
 df=feval('new2',b,z,k); 
 b=b-(f/df); 
 f0=[f0,f]; 
end 
% saídas gráficas 
subplot (1,2,1) 
plot(f0,'m'); 
xlabel('Iteração'); 
ylabel('f(b)'); 
title( ['Flash - Método Newton-Raphson' ] ); 
subplot(1,2,2) 
plot(b0,'m'); 
xlabel('Iteração'); 
63 
 
ylabel('b(% vaporizada da carga)'); 
title( [ 'Flash - Método Newton-Raphson' ] ); 
% salvar dados 
save fbnew.dat f0 -ascii 
save bnew.dat b0 -ascii 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [f]=new1(b,z,k); 
global N 
for i=1:N 
 f(i)=(k(i)-1)*z(i)/(1+b*(k(i)-1)); 
end 
 
f=sum(f); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [df]=new2(b,z,k); 
global N 
for i=1:N 
 df(i)=-(z(i)*(k(i)-1)^2)/(1+b*(k(i)-1))^2; 
end 
df=sum(df); 
11. Aula 11 – Resolução de equação algébrica não linear 
Queridos alunos e alunas, 
Tudo bem? 
Método da Bisseção 
Sistematicamente reduz à metade o intervalo de busca, dado um intervalo inicial 
 de x,x que contém a raiz. 
 
 
 
Recursão do método para resolução de   0xf  
2
xx
x dem

 
Se     0xfxf dm  
Então, me xx  
Caso contrário, md xx  
Crítica do método: 
 fácil de aplicar 
 lento, comparado a outros métodos que veremos a seguir 
65 
 
11.1. Métodos matemáticos para resolução de equações algébricas não lineares 
 
T, P
F, z
V, y
L, x 
 
Flash isotérmico ( psia1600P  e F120T o ) 
Balanço por componente: iii VyLxFz  
Balanço global: VLF  
Restrições: 1z
i
i  , 1x
i
i  , 1y
i
i  
Relação de equilíbrio: iii xky  ,  P,Tkk ii  
Vamos dividir o balanço por componente por F : iii yF
V
x
F
L
z  
Sendo VLF  , 
F
V
1
F
L
 
Vamos definir  como a fração vaporizada da carga (alimentação), então: 
  iii yx1z  
Vamos usar a relação de equilíbrio, então: 
   iii xk1z   
i
i
i k1
z
x

 e 
i
ii
iii k1
zk
xky

 
Vamos usar 1x
i
i  e 1y
i
i  da seguinte forma: 
 
  01k1
z1k
k1
zzk
0xy
i i
ii
i i
iii
i
i
i
i  




 
Portanto a equação de Flash é:      01k1
z1k
f
i i
ii 


 
 
No do 
componente Componente i
z ik 
1 metano 0,7 3,09 
2 propano 0,1 0,39 
3 n-pentano 0,1 0,093 
4 n-hexano 0,1 0,065 
 
11.2. Exemplo: Resolução de equações algébricas não lineares 
11.2.1. Exercício proposto (aula): Programa bis.m e bis1.m 
Implementar o Método da bisseção para o Flash isotérmico; p=1600 psia; T=120°F nas 
condições: 
N Componentes: metano (zM=0,70); propano (zP=0,10); n-pentano (znP=0,10); n-hexano 
(znH=0,10) 
Constante de equilíbrio: metano (kM=3,09); propano (kP=0,39); n-pentano (knP=0,093); n-
hexano (knH=0,065) 
Estimativas iniciais: be=0,20; bd=0,80 
67 
 
 
12. Aula 12 – Resolução de equação algébrica não linear 
 
Queridos alunos e alunas, 
Tudo bem? 
Método da Secante 
Evita o uso de derivadas analíticas, sendo que os pontos empregados para obtenção 
da aproximação numérica da derivada da função sempre cercam a raiz do problema. 
     










1ii
1ii,
numéricai xx
xfxf
xf . 
Da fórmula de Newton-Raphson (  
 i,
i
i1i
xf
xf
xx  ), substituindo  i, xf pela 
aproximação numérica obtemos: 
 
 
 
   












1ii
1ii
i
i,
numéricai
i
i1i
xx
xfxf
xf
x
xf
xf
xx . 
 
f(x)
x
x0
x1Crítica do método da Secante: 
 mais lento que o método de Newton-Raphson 
 porém, evita o uso de derivadas analíticas 
 não há garantia de convergência 
12.1.1. Exemplo de uma função já implementada no OCTAVE 
12.1.1.1. Programa flash1.m 
 
nodo 
componente 
Componente iz ik 
1 metano 0,7 3,09 
2 propano 0,1 0,39 
3 n-pentano 0,2 0,093 
 
% Flash isotérmico; p=1600 psia; T=120°F 
% N Componentes: metano (zM=0,70); propano (zP=0,10); n-pentano (znP=0,20) 
% Constante de equilíbrio: metano (kM=3,09); propano (kP=0,39); n-pentano (knP=0,093) 
% Equação de Flash: f(b)=somatório[(ki-1)zi/(1+b(ki-1))]=0 
% Estimativa inicial: b=0,20 
global N z k f0 
f0=0; 
N=3; 
z(1)=0.70; 
z(2)=0.10; 
z(3)=0.20; 
k(1)=3.09; 
k(2)=0.39; 
k(3)=0.093; 
69 
 
b=0.2; 
bfinal=fzero('flash11',b,1e-4,1); 
f0=f0(2:1:length(f0)) 
feval('flash11',bfinal) 
% saídas gráficas 
plot(f0,'r'); 
xlabel('Iteração'); 
ylabel('f(b)'); 
title(['% vaporizada da carga é ',num2str(bfinal)]), 
% salvar dados 
save fflash.dat f0 -ascii 
save bflash.dat bfinal -ascii 
save fflash f0 
save bflash bfinal 
clear f0 
load fflash.dat 
fflash 
load fflash 
f0 
%************************************************************************** 
function [f]=flash11(b); 
global N z k f0 
for i=1:N 
 f(i)=(k(i)-1)*z(i)/(1+b*(k(i)-1)); 
end 
 
f=sum(f); 
f0=[f0 f]; 
12.1.2. Exercício proposto: Programa scnt.m 
Implementar o Método da Secante para o Flash isotérmico; p=1600 psia; T=120°F, nas 
condições: 
N Componentes: metano (zM=0,70); propano (zP=0,10); n-pentano (znP=0,10); n-hexano 
(znH=0,10) 
Constante de equilíbrio: metano (kM=3,09); propano (kP=0,39); n-pentano (knP=0,093); n-
hexano (knH=0,065) 
71 
 
13. Aula 13 – Resolução de equação algébrica não linear 
Queridos alunos e alunas, 
Tudo bem? 
Método da Regula-Falsi 
Evita o uso de derivadas analíticas, sendo que os pontos empregados para obtenção da 
aproximação numérica da derivada da função sempre cercam a raiz do problema. 
     
ed
ed,
numéricai xx
xfxf
xf


 
Da fórmula de Newton-Raphson (
 
 i,
i
i1i
xf
xf
xx  ), substituindo  i, xf pela aproximação 
numérica obtemos: 
 
 
 
   









ed
ed
i
i,
numéricai
i
i1i
xx
xfxf
xf
x
xf
xf
xx . 
 
f(x)
x
xe
xd
região 
descartada
região 
retida 
 
O critério de descarte de regiões é igual ao critério empregado na Bisseção: 
Se     0xfxf e1i  
Então, 1id xx  
Caso contrário, 1ie xx  
Crítica do método da Regula-Falsi: 
 mais lento que o método da secante 
 porém, evita o uso de derivadas analíticas 
 pode-se garantir a convergência 
13.1.1. Exercício proposto: Programa Regula-Falsi 
Implementar o Método da Regula-Falsi para o Flash isotérmico; p=1600 psia; T=120°F, nas 
condições: 
N Componentes: metano (zM=0,70); propano (zP=0,10); n-pentano (znP=0,10); n-hexano 
(znH=0,10) 
Constante de equilíbrio: metano (kM=3,09); propano (kP=0,39); n-pentano (knP=0,093); n-
hexano (knH=0,065).