Solução do Problema de Pêndulo Oscilante com o Método das Diferenças Finitas

Resumo

 Praticamente todos os engenheiros são confrontados por problemas relacionados ao movimento periódico de corpos livres. Um exemplo simples é um pêndulo oscilante, onde uma partícula de peso W está presa a uma haste sem peso, de comprimento l. Quando a partícula é solta, formando um ângulo inicial θ0 com a vertical, oscila de um lado para o outro com período T [1]. As únicas forças agindo na partícula são a gravidade g e a tensão R na haste, como apresentado na Figura 1. Figura 1: Pêndulo Oscilante Ao aplicar as Leis de Movimento de Newton, deduz-se a seguinte Equação Diferencial Ordinária (EDO) d2θ dt2 g l sin θ  0 , (1) em que θ é o ângulo de deslocamento do pêndulo [1]. Observa-se que a massa m não aparece na equação, pois o movimento de um pêndulo não depende de sua massa. A equação (1) é uma EDO Não-Linear, pela presença do termo sin θ, o que impõe dificuldades para se obter sua solução analítica. Para pequenos deslocamentos angulares, no entanto, sin θ é, aproximadamente, θ, sendo este expresso em radianos. Portanto, nesse caso, a equação (1) se torna linear e assume a forma d2θ dt2 g l θ  0 . (2) bolsista de Iniciação Científica PICME, CNPq/Capes O objetivo deste trabalho é obter as soluções analítica e numérica de um Problema de Valor de Contorno (PVC) para um estudo de caso de um Pêndulo Oscilante. Para isso, considera-se l  0.6096m, g  9.800665m/s2 e as condições de contorno θ(0)  pi/8 e θ(10)  pi/8. Obtém-se então o PVC d2θ dt2  g l θ θ(0)  pi8 e θ(10) pi 8 , (3) cuja solução analítica é θ(t)  1 8 pi ( cos ( 402175607 500 )  1 ) sin ( 402175607x 5000 ) sin ( 402175607 500 )  pi 8 cos ( 402175607x 5000 ) . (4) Para a solução numérica, o algoritmo em MATLAB usado, que implementa o Método das Diferenças Finitas, é dado a seguir. Código 1: Método de Diferenças Finitas - MATLAB syms x %w’’  p(x)w’  q(x)w  r(x), [a,b] dados  inputdlg ({’P: ’, ’Q: ’, ’R: ’, ’n: ’},’Dados ’); limites  inputdlg ({’a: ’, ’f(a): ’, ’b: ’, ’f(b): ’},’PVC’); passo  (str2num(limites {3}) - str2num(limites {1}))/str2num(dados {4}); xVector  str2num(limites {1}) : passo : str2num(limites {3}); %Sistema matricial: Aw  d. Obtendo a matriz tridiagonal A: a(1)  2  (passo 2)*subs(dados {2},x,xVector (2)); %diagonal principal b(1)  -1  (passo /2)*subs(dados {1},x,xVector (2)); %diagonal superior d(1)  -(passo 2)*subs(dados {3},x,xVector (2))  ((1  (passo /2)*subs( dados {1},x,xVector (2)))*str2num(limites {2})); for i  2: length(xVector)-3 a(i)  2  (passo 2)*subs(dados {2},x,xVector(i1)); b(i)  -1  (passo /2)*subs(dados {1},x,xVector(i1)); c(i)  -1 - (passo /2)*subs(dados {1},x,xVector(i1)); %inferior d(i)  -(passo 2)*subs(dados {3},x,xVector(i1)); end a(i1)  2  (passo 2)*subs(dados {2},x,xVector(i2)); c(i1)  -1 - (passo /2)*subs(dados {1},x,xVector(i2)); d(i1)  -(passo 2)*subs(dados {3},x,xVector(i2))  ((1 - (passo /2)*subs (dados{1},x,xVector(i2)))*str2num(limites {4})); %Fatoraçao LU l(1)  a(1); u(1)  b(1)/a(1); z(1)  d(1)/l(1); for i  2: length(xVector)-3 l(i)  a(i) - (c(i)*u(i-1)); u(i)  b(i)/l(i); z(i)  (d(i) - (c(i)*z(i-1)))/l(i); end l(i1)  a(i1) - (c(i1)*u(i)); z(i1)  (d(i1) - (c(i1)*z(i)))/l(i1); %Resolvendo o sistema w(1)  str2num(limites {2}); w(length(xVector))  str2num(limites {4}); w(length(xVector) -1)  z(length(xVector) -2); for i  length(xVector) -2: -1 : 2 w(i)  z(i-1) - (u(i-1)*w(i1)); end A ideia do Método é substituir as derivadas da EDO pelas fórmulas de Diferenças Centradas. O intervalo trabalhado é discretizado, no caso, para [0, 10] foi usado passo h  0.01, formando um sistema de equações de ordem 999  999 para as aproximações θi em cada um dos pontos. O algoritmo gera esse sistema, cuja matriz de coeficientes é tridiagonal e, para sua resolução, aplica Fatoração LU [2]. A Figura 2 apresenta os gráficos de ambas as soluções. Figura 2: Soluções Gráficas Com as soluções gráficas fica claro a proximidade entre os valores numéricos e os analíticos. A Tabela 1 apresenta uma comparação mais clara para alguns dos pontos calculados. Tabela 1: Resultados e Erro i xi θi θ(xi) Erro absoluto 0 0 0.3926990817 0.3926990817 0 200 2 0.9313500208 0.9445509565 0.0132009357 400 4 -0.7001842363 -0.7081023307 0.0079180944 600 6 -0.7001842363 -0.7081023307 0.0079180944 800 8 0.9313500208 0.9445509565 0.0132009357 1000 10 0.3926990817 0.3926990817 0 O método provou sua eficiência, pois realmente aproximou a solução numérica à analítica.