Utilizar as funções nativas do MATLAB® é sempre a melhor opção?
“Everything should be made as simple as possible, but no simpler.” – Einstein
Você consegue se enxergar cenário a seguir?
Digita-se tudo corretamente, segue-se a lógica adequada, clica-se no “run” e... “ufa, o código rodou”.
Todos nós, estudantes de engenharia e engenheiros, passamos por este pequeno momento glorioso durante nossas soluções de problemas utilizando as mais diversas linguagens de programação, em especial o MATLAB™. Cá entre nós, depois de uma trabalheira danada, nada melhor que ver nosso código rodando.
Porém, logo após esta felicidade fugaz, vem outro desafio: será que a solução está correta? O gráfico, o vetor ou a matriz está condizente com o esperado? A resposta faz sentido? Em caso negativo, alteramos o método, implementamos algumas melhorias, o código roda, porém mais uma vez a resposta continua não fazendo sentido. Temos certeza que fizemos tudo certo, então, talvez nossa inclinação seja simplesmente desistir. Afinal, se é tão difícil para o MATLAB™, nenhuma pessoa razoável pode esperar que obtivéssemos uma solução. Porém, como bons futuros engenheiros, não procuramos desculpas e nosso único recurso é desenvolver uma solução com base em nossos conhecimentos de métodos numéricos.
Antes de tudo, devemos ter em mente que nunca podemos confiar cegamente em todo resultado gerado por um computador. Existe uma velha máxima que diz: “se entrar lixo, sai lixo”. Infelizmente, algumas pessoas pensam que, independentemente do que entrou (os dados) e do que está acontecendo no interior (o algoritmo), o computador gera um resultado religiosamente correto, o que não é verdade.
Para que evitemos cair neste tipo de armadilha, vejamos um exemplo a seguir de como resolver numericamente uma EDO, mesmo que inicialmente o problema pareça não ter solução:
Enunciado. Uma vazão de água entra em um
taque cilíndrico a uma taxa constante Qentrada, e enche o tanque até
atingir o nível yalto. Neste ponto, a água escoa para fora do tanque
através de um tubo de descarga circular (um sifão), produzindo uma fonte na
saída do tubo. A fonte funciona até que o nível de água diminua para ybaixo;
depois disso o sifão enche-se com ar e a fonte para. O ciclo é então repetido à
medida que o tanque enche-se de água até atingir yalto e a fonte
jorra de novo.
Quando o sifão está funcionando,
o fluxo para fora Qsaída pode ser calculado com a fórmula baseada na
Lei de Torricelli: Qsaída = C*√(2*g*y)*ℼ*r^2 (1)
 |
| FIGURA 1 - Tanque cilíndrico com sifão |
Desprezando o volume de água no
tubo, calcule e trace o nível de água no tanque como uma função do tempo durante
100s. Considere uma condição inicial de um tanque vazio y(0) = 0 e empregue os
parâmetros seguintes para seus cálculos:
Rt = 0,05m; yalto
= 0,1m; Qentrada = 50x10^-6 m^3/s; r = 0,007 m; C = 0,6; ybaixo
= 0,025m; g = 9,81 m/s^2
Solução:
Quando a fonte está funcionando,
a taxa de variação no volume V (m^3) do tanque é determinada por um balanço
simples do fluxo que entra menos o fluxo que sai do tanque:
dV/dt = Qentrada - Qsaída (2)
Como o tanque é cilíndrico, V = ℼ*Rt^2*y. Substituindo essa relação
em (1) e em (2), obtemos:
dy/dt = [Qentrada -C*√(2*g*y)*ℼ*r^2]/ℼ*Rt^2*y
Quando a fonte não está
funcionando, o segundo termo no numerador vai à zero. Podemos incorporar esse
mecanismo no modelo introduzindo uma nova variável adimensional nomeada sifão, que é igual a zero, quando a
fonte está parada, e igual a um, quando a fonte está funcionando:
dy/dt = [Qentrada - sifao* C*√(2*g*y)*ℼ*r^2]/ℼ*Rt^2*y
Neste contexto, a
variável sifão pode ser pensada como
um interruptor que “liga” e “desliga” a fonte. Essas variáveis de dois estados
são chamadas de variáveis lógicas ou
booleanas, nas quais zero é
equivalente a falso, e um é equivalente a verdadeiro.
Em seguida, devemos
relacionar sifão à variável
dependente y. A princípio, a variável sifão
é definida como zero sempre que o nível cair abaixo de ybaixo e como
um sempre que o nível subir acima de yalto.
A função a seguir do
MATLAB™ segue essa lógica no cálculo da derivada:
function dy =
ex_edo(t,y)
global sifao
Rt = 0.05; r = 0.007; yalto =
0.1; ybaixo = 0.025; C = 0.6; g = 9.81; Qentrada = 0.00005;
if y(1) <= ybaixo
sifao = 0;
elseif y(1) >=
yalto
sifao = 1;
end
Qsaida =
sifao*C*sqrt(2*g*y(1))*pi*r^2;
dy = (Qentrada-Qsaida)/(pi*Rt^2);
Observe que, como seu
valor deve ser mantido entre chamadas da função, a variável sifao é declarada como uma variável
global. Embora o uso de variáveis globais não seja indicado (particularmente em
programas extensos), esse uso é útil no contexto presente.
O programa a seguir
emprega a função nativa ode45* para integrar ex_edo e gerar um gráfico da
solução:
clc; close all; clear all;
global sifao
sifao = 0;
tspan = [0 100]; y0 = 0;
[tp,yp] =
ode45(@ex_edo,tspan,y0);
plot(tp,yp);
xlabel('Tempo (s)')
ylabel('Nivel de agua
no tanque (m)')
* ode45 – A função
ode45 usa um algoritmo que usa simultaneamente fórmulas de Runge-Kutta de
quarta e quinta ordens para resolver a EDO e fazer estimativas do erro para
ajuste do tamanho do passo. O MATLAB™ recomenda que a ode45 é a melhor função
para aplicar como uma “primeira tentativa” para a maioria dos problemas.
Como mostra a Figura 2,
o resultado está claramente incorreto. Exceto para o processo de enchimento
original, o nível parece começar a decrescer antes de atingir yalto.
De modo similar, quando o sifão está jorrando água, ele é fechado bem antes do
nível cair para ybaixo.
 |
| FIGURA 2 - Plot do código com função nativa do MATLAB |
Neste ponto,
suspeitando que o problema demande uma função computacional mais poderosa que a
rotina confiável ode45, você pode ser tentado a usar algum dos outros
solucionadores de EDOs do MATLAB, tais como ode23s ou ode23tb. Porém se você
fizer isso, descobrirá, que embora rotinas produzam resultados um pouco
diferentes, elas ainda geram solução incorretas.
A dificuldade surge
porque a EDO é descontinua no ponto em que o sifão é desligado ou ligado. Por
exemplo, quando o tanque está enchendo, a derivada depende apenas do fluxo que
entra no tanque o e, para os parâmetros presentes, tem um valor constante de
6,366*10^-3 m/s. Entretanto, assim que o nível atinge yalto, o fluxo
que sai entra em ação e a derivada cai de forma abrupta para -1,013*10^-2 m/s.
Embora rotinas com controle adaptativo do tamanho do passo usadas pelo MATLAB™
funcionem de modo eficaz para vários problemas, elas muitas vezes têm problemas
quando lidam com essas descontinuidades.
Como o problema resulta
da aplicação de um tamanho de passo adaptativo por meio de uma descontinuidade,
pode-se reverter o método utilizando um tamanho de passo pequeno e constante
(que é uma abordagem mais simples). Podemos implementar essa estratégia de
solução substituindo ode45 por uma outra função, de passo constante, criada
manualmente. Vamos nomeá-la de rk4sys.
function [tp,yp] =
rk4sys(dydt,tspan,y0,h,varargin)
% rk4sys: Runge-Kutta de
quarta ordem para um Sistema de EDOs
% [t,y] =
rk4sys(dydt,tspan,y0,h,p1,p2,...): integra um sistema de EDOs com o método RK
de quarta ordem
% entrada:
% dydt = nome do arquivo-M
que avalia as EDOs
% tspan = [ti,tf];
instantes inicial e final com a saída geradas em intervalos de h, ou = [t0 t1
...tf]; em instantes específicos
% h = tamanho do passo
% p1,p2,... = parâmetros
adicionais usados por dydt
% saída:
% t = vetor da variável
independente
% y = vetor de solução para
a variável dependente
if nargin<4, error('são
necessários pelo menos 4 argumentos de entrada')
end
if any(diff(tspan)<=0), error('tspan não está
em ordem crescente'),
end
n = length(tspan);
ti = tspan(1); tf = tspan(n);
if n == 2
t = (ti:h:tf)'; n = length(t);
if t(n)<tf
t(n+1) = tf;
n = n+1;
end
else
t = tspan;
end
tt = ti; y(1,:) = y0;
np = 1; tp(np) = tt; yp(np,:) =
y(1,:);
i = 1;
while(1)
tend = t(np+1);
hh = t(np+1) - t(np);
if hh>h,hh = h;
end
while(1)
if tt+hh>tend,hh = tend-tt;
end
k1 = dydt(tt,y(i,:),varargin{:})';
ymid = y(i,:) + k1.*hh./2;
k2 = dydt(tt+hh/2,ymid,varargin{:})';
ymid = y(i,:) + k2.*hh./2;
k3 = dydt(tt+hh/2,ymid,varargin{:})';
yend = y(i,:) + k3.*hh;
k4 = dydt(tt+hh,yend,varargin{:})';
phi = (k1+2*(k2+k3)+k4)/6;
y(i+1,:) = y(i,:) + phi*hh;
tt = tt+hh;
i=i+1;
if tt>=tend,break,end
end
np = np+1; tp(np) = tt; yp(np,:) = y(i,:);
if tt>=tf,break,end
end
Para o programa rodar
só é necessário substituir a quarta linha do código apresentado anteriormente por:
clc; close all; clear all;
global sifao
sifao = 0;
tspan = [0 100]; y0 = 0;
[tp,yp] =
rk4sys(@ex_edo,tspan,y0,0.0625);
plot(tp,yp);
xlabel('Tempo (s)')
ylabel('Nivel de agua
no tanque (m)')
O resultado obtido será:
 |
| FIGURA 3 - Plot do código com função criada rk4sys |
Agora sim a solução evolui como esperado. O tanque enche até yalto e depois
esvazia até atingir ybaixo,quando o ciclo se repete.
Situações como aquela
descrita na Figura 1 são especialmente perigosas – isto é, embora a saída
seja incorreta, ela não é obviamente incorreta; ou seja, a simulação não se
torna instável ou produz níveis negativos. Na verdade, a solução se move para
cima e para baixo na forma de uma fonte intermitente, embora incorreta.
Esperamos que este
estudo de caso ilustre que mesmo uma grande parte de pacotes computacionais
como o MATLAB™ não é infalível e nem sempre trazem funções nativas que
necessitamos. Portanto, estudantes de engenharia e engenheiros sempre examinam
saídas numéricas com um ceticismo saudável com base em suas experiências e
conhecimentos dos problemas que estão resolvendo.
_________________________________________________________________________________
Este problema foi baseado no estudo de caso contido em:
CHAPRA, Steven C. Métodos numéricos aplicados com MATLAB para engenheiros e cientistas; tradução Rafael Silva Alípio. - 3. ed. - Porto Alegre: AMGH, 2013. Pgs 609-612.
O livro a seguir foi de grande ajuda também:
GILAT, Amos. MATLAB com aplicações em engenharia; tradução Rafael Silva Alípio. 4. ed. - Porto Alegre: Bookman, 2012.
Caso tenha ficado alguma dúvida ou tenha encontrado algum erro, sinta-se a vontade para enviar um e-mail para nós colaboradores do blog.
Filipe Natã
Comentários
Postar um comentário