Exercícios de Cálculo Numérico em Scilab/Matlab

Prévia do material em texto

Roteiro de exercícios - Aula prática 4
Instruções Gerais
A entrega dos exercícios pode ser feita de forma individual ou em dupla (OU SEJA, GRUPOS
COM 3 OU MAIS ALUNOS NÃO SERÃO ACEITOS). Cada grupo deverá entregar um único
relatório contendo suas respostas às perguntas dadas a seguir, bem como os arquivos com os códigos
desenvolvidos. Entregas do relatório sem os códigos ou dos códigos sem o relatório receberão nota
0. Não é necessário fazer um relatório completo, com introdução, metodologia e afins, basta
responder cada um dos itens. Os arquivos devem ser zipados e enviados via Moodle até as 23h59
do dia 08/12.
Os arquivos (códigos e relatório) serão testados por software detector de plágio. Caso se con-
firme a cópia, será atribuída nota 0 a todos os grupos envolvidos (O argumento �eu fiz o trabalho
e emprestei para eles darem uma olhada, não sabia que eles iam copiar� não é aceitável).
Neste laboratório serão utilizados dados relativos a um experimento feito com um reator tubular
(PFR, plug flow reactor). Diz-se que o escoamento nesse reator é do tipo pistonado, vide por
exemplo https://www.youtube.com/watch?v=9LkSl93nKKA.
Os dados de um experimento análogo ao do vídeo estão disponíveis na planilha PFR.xls que
está no Moodle. Nesse experimento, uma injeção de um composto traçador é feita na entrada do
reator e sua concentração na saída é medida. A partir dos dados do arquivo PFR.xls, pode-se gerar
o seguinte gráfico da concentração do traçador na saída do reator em função do tempo:
O traçador não reage com o solvente, apenas escoa pela tubulação. No início do experimento,
a concentração na saída é 0 porque o traçador �ainda não chegou� até o final. Após um certo
tempo, o traçador começa a sair mas, devido à difusão interna dentro do reator, ele não vai sair
todo de uma vez. Repare no vídeo anterior que o setor vermelho corresponde ao ponto com maior
concentração do traçador e que há um gradiente de concentração do traçador dentro do reator. Em
termos técnicos,diz-se que há uma distribuição do tempo de residência, ou seja, há uma distribuição
de probabilidade em que é improvável que haja uma concentração alta de traçador na saída para
instantes de tempo muito pequenos ou muito grandes.
Para ler os dados do Excel, utilize os seguintes comandos para salvar os dados de tempo (s)
contidos na primeira coluna em um vetor chamado t e os dados de concentração (mol/L) da segunda
coluna em um vetor chamado c. O comando plot serve apenas para testar que a importação ocorreu
com sucesso.
sheets = readxls('PFR.xls')
s1 = sheets(1);
t = s1.value(2:\$,1); // excluindo a primeira linha porque é o cabeçalho
c = s1.value(2:\$,2); // excluindo a primeira linha porque é o cabeçalho
plot(t,c,'.')
xlabel('Tempo(s)')
ylabel('Concentração (mol/L)')
Se utilizar Matlab, os comandos apropriados são
[Num,txt]=xlsread('PFR.xls') % separando a parte numérica da parte de texto
t = Num(:,1);
c = Num(:,2);
plot(t,c,'.')
xlabel('Tempo(s)')
ylabel('Concentração (mol/L)')
Nota: é preciso que o arquivo do Scilab/Matlab esteja na mesma pasta que o arquivo .xls. Além
disso, é preciso que a pasta de operação do programa seja a mesma que contém esses arquivos:
O tempo de residência de um reator é definido como
τ =
∫∞
0
tc(t)dt∫∞
0
c(t)dt
em que t é o tempo e c é a concentração. Nos próximos exercícios, admite-se que a concentração
nos últimos pontos experimentais já é suficientemente próxima de 0, de modo que o limite superior
da integral pode ser tomado como o final do experimento, ao invés de considerar um tempo infinito.
1. A função inttrap(x,y) do Scilab (trapz(x,y) no Matlab) calcula a integral por trapézios
a partir dos vetores x e y, que correspondem à variável independente e à função nesses pontos,
respectivamente. Calcule o valor aproximado do tempo de residência desse reator utilizando o
método dos trapézios.
2. Crie uma função chamada intSimpson com o seguinte cabeçalho,
function integral = intSimpson(h,y)
O vetor y corresponde a f(x) em pontos igualmente espaçados no eixo x com distância h entre
eles. O retorno da função, a variável chamada integral, é o valor aproximado da integral de
acordo com a regra de Simpson de 1/3:
∫ b
a
f(x)dx ≈ h
3
f(x0) + 4
n−1∑
i = 1
i ímpar
f(xi) + 2
n−2∑
i = 2
i par
f(xi) + f(xn)

Utilize a sua função para calcular uma aproximação do tempo de residência pelo métodod e
Simpson. Compare com o resultado do item anterior. Dicas: tome cuidado com o fato de que a
numeração de vetores em Scilab/Matlab começa em 1. Para calcular o módulo da divisão de um
número n por um número q, utilize modulo(n,q) em Scilab e mod(n,q) em Matlab.
3. Crie uma função com o seguinte cabeçalho,
function dydx = derivada(y, i, flag)
y é um vetor de observações da variável dependente, i é a posição em que se deseja calcular a
derivada e flag indica qual tipo de derivada que se deseja calcular (flag = 1 corresponde a deri-
vada para trás, flag = 2 corresponde a derivada para a frente e flag = 3 corresponde a derivada
centrada). Utilize aproximações que envolvam apenas dois pontos do vetor y a cada vez. Preveja
em seu programa os casos em que o input do usuário pode causar problemas ao tentar acessar
posições fora do vetor. Avalie os três tipos de derivada no instante t = 170 (correspondente a i =
86).
4. Um método alternativo para o cálculo de derivadas consiste em encontrar por mínimos quadra-
dos um polinômio que se ajusta aos dados e então calcular a derivada desse polinômio no ponto
desejado. Escolha uma vizinhança apropriada em torno do ponto t = 170 (correspondente a i =
86) e estime a derivada nesse ponto usando como polinômio aproximador uma reta e uma parábola
(um caso de cada vez). Utilize a função para mínimos quadrados desenvolvida no Laboratório 3.
Para facilitar o cálculo da derivada, crie uma função com o seguinte cabeçalho
function dydx = derivada(c, x)
Em que c é um vetor com os coeficientes do polinômio e x é um escalar correspondente ao ponto
em que se deseja calcular a derivada. Faça sua função de forma genérica, permitindo que c possa
ter qualquer tamanho. Dica: tome cuidado em relação à ordem dos coeficientes de c, ou seja, c(1)
corresponde ao coeficiente de qual potência de x?
5. Compare os resultados obtidos nas questões 3 e 4. Qual parece ser mais adequado para o
caso em questão?