Análise de Sistemas Stiff com Método de Euler e Funções MATLAB

Sistemas Stiff e descontinuos Suponha que aplicamos o método de Euler para resolver o mesmo problema numericamente a O método de Euler explícito resulta b O método de Euler implícito resulta Exemplo Considere agora um sistema linear de EDO dy1dt 5y1 3y2 dy2dt 100y1 301y2 Aplicando o método de Euler implícito y1i1 y1i 5y1i1 3y2i1h y2i1 y2i 100y1i1 301y2i1h que pode ser reescrito 1 5hy1i1 3y2i1 y1i 100y1i1 1 301hy2i1 y2i Notase que para encontrar a solução é preciso resolver um conjunto de equações em cada passo i Se o sistema de EDO para não linear a solução se torna ainda mais complexa pois envolve resolver um sistema de equações simultâneas não lineares Assim embora a estabilidade da solução seja obtida por meio de abordagens implícitas há um preço a ser pago pelo aumento da complexidade da solução Funções MATLAB para sistemas stiff O MATLAB possui quatro funções ode15s ode23s ode23t ode23tb para resolver sistemas rígidos de ODEs enquanto que o Octave apenas a ode15s Exemplo A equação de Balthasar van der Pol é um modelo de circuito eletrônico da época dos tubos de vácuo d²y2dt² μ1 y1²dy1dt y1 0 A solução para esta equação se torna progressivamente mais rígida à medida que aumenta Dadas as condições iniciais y10 1 e y10 1 resolve os dois casos a seguir a para μ 1 use a função ode45 para encontrar a solução desde t 0 até 20 e b para μ 1000 use ode15s para resolver desde t 0 até 6000 A primeira tarefa é representar a equação de van der Pol no espaço de estados dy1dt y2 dy2dt μ1 y1²y2 y1 As funções vanderpolfunmu1m e vanderpolfunmu1000m foram criadas function yp vanderpolfunmu1ty mu 1 yp y2mu1y12y2y1 function yp vanderpolfunmu1000ty mu 1000 yp y2mu1y12y2y1 A chamada da função ode45 Para μ 1 o sistema não é stiff ty ode45 vanderpolfunmu10 201 1 subplot121 plotty1 subplot122 plotty2 A chamada da função ode15s Para μ 1000 o sistema é stiff ty ode15s vanderpolfunmu10000 60001 1 subplot121 plotty1 subplot122 plotty2 Descrição das Funções MATLAB para sistemas stiff ode15s Esta função é um solver de ordem variável que utiliza fórmulas de diferenciação numérica É um solver de passo múltiplo que opcionalmente usa as fórmulas de diferenciação reversa Ela é recomendada para problemas rígidos de baixa a média precisão também disponível no Octave ode23s Esta função é baseada em uma fórmula modificada de Rosenbrock de ordem 2 Por ser um solver de passo único pode ser mais eficiente que o ode15s Ela pode resolver alguns tipos de problemas rígidos melhor que os ode15s ode23t Esta função implementa a regra trapezoidal com um interpolante livre Ela é usada para problemas moderadamente rígidos com baixa precisão onde você precisa de uma solução sem amortecimento numérico ode23tb Esta é uma implementação de uma fórmula implícita de RungeKutta com um primeiro estágio que é uma regra trapezoidal e um segundo estágio que é uma fórmula de diferenciação reversa da ordem 2 Este solver também pode ser mais eficiente que o ode15s Seleção do solver de ODE e DAE do MATLAB httpswwwmathworkscomhelpmatlabmathchooseanodesolverhtml httpsblogsmathworkscomloren20150923odesolverselectioninmatlabstidsrchtitle Accuracy Problemas descontínuos O filósofo natural Romano Plínio o Velho supostamente tinha um chafariz intermitente em seu jardim como o da figura A água entra em um tanque cilíndrico com uma vazão constante Qin e enche até que a água alcance o nível yhigh Nesse ponto a água sai do tanque através de um tubo de descarga circular produzindo uma fonte na saída do tubo O chafariz funciona até que o nível da água diminua até ylow quando o sifão se enche de ar e a fonte para O ciclo então se repete sucessivamente Qout cπr²2gy Quando o sifão está em ação a saída Qout pode ser calculada pela lei de Torricelli Qout cπr²2gy Desprezando o volume de água no tubo calcule e faça o gráfico do nível da água no tanque em função do tempo por 100 segundos Suponha que inicialmente o tanque está vazio y0 0 e utilize os seguintes dados RT 005 m r 0007 m ylow 0025 m yhigh 01 m C 06 g 981 ms² Qin 50 x 10⁶ m³s Aplicando o balanço volumétrico na fase em que há fluxo através do sifão ou seja para y yhigh dydt Qin Qout dydt Qin cπr²2gyπR² Quando não há fluxo através do sifão y ylow o segundo termo não está presente Para usarmos a mesma EDO para as duas fases podese acrescentar uma variável lógica binária siphon dydt Qin siphoncπr²2gyπR² No código da ODE devemos relacionar o valor da variável lógica siphon à variável dependente y Observe que como seu valor deve ser mantido entre as chamadas da função a siphon foi declarada como variável global function dy Pliniodefunt y global siphon Rt 005 r 0007 yhi 01 ylo 0025 C 06 g 981 Qin 000005 if y1 ylo siphon 0 elseif y1 yhi siphon 1 end Qout siphon C sqrt2 g y1 pi r 2 dy Qin Qout pi Rt 2 Chamada da função ode45 no Octave Plinyoderun global siphon siphon 0 tspan 0 100 y0 0 tpyp ode45 Plinyodefuntspany0 plottpyp xlabeltime s ylabelwater level in tank m Note pelo gráfico que resultado está claramente incorreto Tentando também a função ode15s no Octave obtemos solução também incorreta A dificuldade surge porque a ODE é descontínua no ponto em que o sifão liga ou desliga Embora as rotinas de passo adaptativo do MATLAB funcionem bem para muitos problemas elas podem não lidar bem com descontinuidades A solução é usar um passo pequeno e constante Podemos implementar essa estratégia de solução simplesmente substituindo o ode45 pela função rk4fun e rk4run que nós criamos trocando a função dydtfun pela Plinyodefun e usando o passo h 00625 que é suficientemente pequeno Exemplo Considere o problema completo do bungee jumper e determine a posição e a velocidade do jumper para os seguintes parâmetros L 30 m g 981 ms² m 681 kg cd 025 kgm k 40 Nm e γ 8 Nsm Encontre as soluções desde t 0 até 50 s e assuma que as condições iniciais são x0 v0 0 A equação do movimento dvdt g cdm sgnvv² Quando o cabo começa a esticar x L a lei de Hooke pode ser usada para representar o efeito elástico do cabo Além disso uma força de amortecimento também deve ser incluída para levar em consideração os efeitos de atrito à medida que o cabo se estica e se contrai Nessa fase a EDO muda para incluir novas forças resistentes dvdt g cdm sgnvv² kmx L γmv Como é necessário saber a posição do jumper incluiuse a equação diferencial para a distância dxdt v A função bungeefunm function dydt bungeefuntyLcdmkgamma g 981 cord 0 if y1 L se o cabo exerce força cord kmy1Lgammamy2 end dydt y2 g signy2cdmy22cord 238 The following nonlinear parasitic ODE was suggested by Hornbeck 1975 dydt 5y t2 If the initial condition is y0 008 obtain a solution from t 0 to 5 a Analytically b Using the fourthorder RK method with a constant step size of 003125 c Using the MATLAB function ode45 d Using the MATLAB function ode23s e Using the MATLAB function ode23tb Present your results in graphical form