Uma exposição didática do método de Runge-Kutta no estudo do pêndulo amortecido

A didactic exposition of the Runge-Kutta method in the study of the damped pendulum

Resumos

Diversos fenômenos naturais são frequentemente modelados usando equações diferenciais. Contudo, uma solução analítica de tais equações nem sempre é possível, de forma que métodos numéricos devem fazer parte do arsenal dos físicos para estudar seu sistema de interesse. Nesse trabalho esperamos proporcionar uma apresentação pedagógica para projetos que necessitem encontrar as soluções de equações diferenciais ordinárias por métodos numéricos. Para isso, usamos o método de Runge-Kutta de quarta ordem para descrever o pêndulo simples com a presença de uma força resistiva através dos seus diagramas de fase nos amortecimentos subcrítico, crítico e supercrítico, por fim comparamos as soluções numéricas com os resultados analíticos.

Palavras-chave:
Métodos numéricos; Runge-Kutta; Pêndulo amortecido

Various natural phenomena are often modeled using differential equations. However, an analytical solution of such equations is not always possible, so numerical methods must be part of the physicists’ arsenal to study their system of interest. In this work we hope to provide a pedagogical presentation for projects that need to find solutions to ordinary differential equations using numerical methods. For this, we use the fourth-order Runge-Kutta method to describe the simple pendulum with the presence of a resistive force through its phase diagrams in subcritical, critical and supercritical damping, finally we compare the numerical solutions with the analytical results.

Keywords:
Numerical methods; Runge-Kutta; Damped pendulum

1. Introdução

Equações diferenciais têm um papel extremamente importante em todas as áreas da Física. Contudo, uma solução analítica pode ser extremamente difícil, seja pela dificuldade matemática, fugindo do escopo do projeto, seja pelo tempo ou aproximações necessárias para sua resolução. Existem ainda os casos em que soluções analíticas são impossíveis de serem obtidas. Nesse cenário, os métodos numéricos surgem como uma ferramenta poderosa nas mãos de pesquisadores e estudantes. Existem diversos métodos numéricos, como o método de Euler, Metropolis, Decomposição de Adomian e muitos métodos especiais são desenvolvidos ainda hoje no estudo de casos específicos. Métodos numéricos computacionais tomaram grandes proporções com a popularização dos computadores portáteis. A escolha de um método depende do problema abordado e do equilíbrio que o pesquisador deseja entre o tempo de computação e a precisão do algoritmo. Dentre eles, um dos mais conhecidos são os métodos de Runge-Kutta, que foram desenvolvidos por volta de 1895 pelos matemáticos alemães Carl David Tolmé Runge e Martin Wilhelm Kutta. O método tradicionalmente é utilizado pela sua eficácia em resolver equações diferenciais de primeira ordem com problema de valor inicial. Consiste em uma família de métodos iterativos para aproximar numericamente a solução exata de uma equação diferencial ordinária pelos primeiros n termos da expansão em série de Taylor, aproveitando o melhor do método de Euler [¹].

No presente trabalho, focaremos no método de Runge-Kutta com passos de tamanho fixo. Esse método possui certos inconvenientes, como o gasto computacional em regiões da função estudada onde passos maiores seriam suficientes para obter soluções com erros pequenos. Apesar de existirem diversos métodos que adaptam o tamanho do passo, pretendemos apresentar uma introdução ao método, por isso, estudaremos o método mais simples, com tamanho de passo fixo.

Utilizamos como sistema para aplicar o método numérico de Runge-Kutta de quarta ordem o pêndulo simples com amortecimento. Há mais de 400 anos o pêndulo é um dos instrumentos mais estudados na física, desde as observações de Galileu nos séculos XVI e XVII, fundamentais para o desenvolvimento de sua cinemática e da dinâmica newtoniana [²]. Equações de movimento [³,⁴,⁵], efeitos de memória e não-linearidades [⁶, ⁷], correções na amplitude [⁸] e período [⁹], e grandezas medidas em regimes de grandes oscilações [¹⁰] continuam como temas de grande interesse. Por esta razão, a descrição de seu movimento está presente na maioria dos livros-texto de física, seja de cursos superiores ou do ensino médio. Além do interesse no próprio sistema, sua relevância é ampliada por muitos fenômenos físicos serem governados por equações diferenciais semelhantes, como processos de relaxação [¹¹], elásticas [¹²], junções Josephson bosônicas [¹³], entre outros.

Assim, nosso objetivo é apresentar uma forma sistemática de implementar o método de Runge-Kutta de quarta ordem de forma simples e didática. Mais ainda, utilizamos o método para obter o diagrama de fases do pêndulo em casos onde existe a presença de uma força resistiva. Para isso, o texto a seguir é dividido em quatro seções. Na segunda seção, descrevemos o método de Runge-Kutta de quarta ordem e sua utilização na solução de uma equação diferencial. Na terceira seção, discutimos o sistemas físico no qual o método foi usado para obter os seus diagramas de fase. Nossos resultados numéricos são apresentados na quarta seção e os comparamos com as soluções analíticas. Por fim, sumarizamos os resultados e apresentamos as considerações finais nas conclusões.

2. Runge-Kutta de Quarta Ordem

Os métodos de Runge-Kutta foram desenvolvidos por volta de 1895 e possuem grande eficácia em resolver equações diferenciais de primeira ordem com problema de valor inicial. O objetivo dessa seção é demonstrar as equações do método e sua implementação em um algoritmo computacional. Iniciaremos descrevendo o método de Euler, que corresponde ao método de Runge-Kutta de primeira ordem. Para isso, considere a equação diferencial linear de primeira ordem não-homogênea

$$\frac{dy(x)}{dx}+k(x)y(x)=w(x),$$

que de maneira geral podemos também escrever como

(1)$$y^{\prime }(x)=F(x,y),$$

com y′(x) = dy(x)/dx e o valor inicial dado por y(x₀) = y₀. Queremos encontrar uma aproximação para a solução y(x) dessa equação diferencial a partir do valor inicial. Para isso, começamos considerando a expansão em série de Taylor de uma dada função y = f(x) próxima a um ponto a, dada por [¹⁴],

(2)$$y(x)=y(a)+(x-a)y^{\prime }(a)+\frac{{(x-a)}^2}{2!}{y}^{″}(a)+\frac{{(x-a)}^3}{3!}{y}^{‴}(a)+\cdots ={\displaystyle \sum _{m=0}^{\infty }\frac{{(x-a)}^m}{m!}{y}^{(m)}(a)}.$$

A expansão em série de Taylor nos permite avaliar funções na proximidade de determinado ponto, como o ponto a acima. Assim, na proximidade de x₀, o ponto x é escrito como x = x₀ + h, onde h é chamado de tamanho do passo, ou seja, o ponto x onde estamos avaliando a função está a uma distância h do ponto inicial. Com isso, a partir da equação (2), a expansão em torno de x₀ pode ser escrita na forma

(3)$$y(x)=y(x_0+h)=y(x_0)+h{y}^{\prime }(x_0)+\frac{h^2}{2!}{y}^{″}(x_0)+\frac{h^3}{3!}{y}^{‴}(x_0)+\cdots$$

Na expansão acima vamos considerar o tamanho do passo pequeno, de tal forma que podemos desprezar os termos de ordens superiores (h₂, h₃, ...). Portanto, o método de Euler consiste em aproximar a função y(x) da seguinte forma

(4)$$y(x_0+h)\simeq y(x_0)+h{y}^{\prime }(x_0).$$

2.1. Implementação

2.1.1. Métodos de Euler e de Runge–Kutta de segunda ordem

O domínio da função, o eixo x, deve ser dividido em n partes ($n\in \mathbb{N}$). Toma-se um ponto inicial x₀ e calcula a função no ponto subsequente, ou seja, x₁ = x₀ + 1 × h usando a equação (4). A partir de x₁ calculamos a função em x₂ = x₁ + h = x₀ + 2 × h. Portanto, o valor para a função num ponto posterior x_n+1 = x_n + h é obtido a partir do ponto anterior x_n = x₀ + nh usando a aproximação (com erro da ordem h²):

(5)$$y(x_{n+1})\simeq y(x_n)+h{y}^{\prime }(x_n),$$

que traduz o método de Euler ou método da reta tangente. Com o propósito de verificar a precisão dessa aproximação consideremos como exemplo a função exponencial f(x) = e^x com o tamanho do passo h = 0, 1. A tabela 1 mostra a comparação dos valores reais desta função em três valores de x₀ (0, 5, 1, 0 e 2, 0) com a aproximação pelo método de Euler. O erro absoluto, $‖\text{valor real - valor aproximado}‖$, também foi calculado e mostrado na tabela 1 para n = 1 (um passo) e n = 10 (dez passos).

Tabela 1
Valores exatos e aproximados pelo método de Euler da função e^x e os erros absolutos em três valores de x₀ com h para n = 1 e n = 10.

x 0	xn = x0 + nh (n = 1)	valor exato f(x1)	valor aproximado f(x1)	erro absoluto
0,5	0,6	1,82212	1,81359	0,00853
1,0	1,1	3,00417	2,99011	0,01406
2,0	2,1	8,16617	8,12796	0,03821
x 0	xn = x0 + nh (n = 1)	valor exato f(x1)	valor aproximado f(x1)	erro absoluto
0,5	1,5	4,48169	4,46072	0,02097
1,0	2,0	7,38906	7,35448	0,03457
2,0	3,0	20,08554	19,9916	0,09398

Vemos que o erro cresce quando aumentamos o número de passos, ou seja, o número de passos entre a condição inicial e o ponto em que estamos calculando a função. Apesar do erro do truncamento diminuir quando usamos um valor h menor, o erro do arredondamento acaba aumentando. Para reduzir o erro, faremos uma nova aproximação.

Considere agora, usando a equação (2), expansões em série de Taylor da função y(x) em torno do ponto x_n+1/2 = x_n + h/2, com um tamanho de passo de h/2. A primeira expansão é aplicada a y(x_n+1), produzindo

(6)$$y(x_{n+1})=y(x_n+h/2)+\frac{h}{2}{y}^{\prime }(x_n+h/2)+\frac{1}{2!}{\left(\frac{h}{2}\right)}^2{y}^{″}(x_n+h/2)+\mathcal{O}(h^3),$$

em que $\mathcal{O}\text{(}{h}^{\text{3}}\text{)}$ indica termos envolvendo potências de h com expoente maior ou igual a 3. A segunda expansão é aplicada “para trás”, a y(x_n), produzindo

(7)$$y(x_n)=y(x_n+h/2)-\frac{h}{2}{y}^{\prime }(x_n+h/2)+\frac{1}{2!}{\left(\frac{h}{2}\right)}^2{y}^{″}(x_n+h/2)+\mathcal{O}(h^3),$$

Subtraindo a equação (7) da equação (6), obtemos

(8)$$y(x_{n+1})-y(x_n)=h{y}^{\prime }(x_n+h/2)+\mathcal{O}(h^3),$$

já que os termos proporcionais a h² se cancelam. Com isso, obtemos uma aproximação correta até ordem h³:

(9)$$y(x_{n+1})\simeq y(x_n)+h{y}^{\prime }(x_n+h/2).$$

A equação acima é conhecida como método de Euler modificado, método de Heun ou método de Runge-Kutta de segunda ordem. Comparando os resultados das expressões (5) e (9), observamos que a última aproximação apresenta um resultado mais preciso, uma vez que o valor da derivada é calculada num ponto mais próximo de x₀. O erro envolvido nesse caso é da ordem h³.

2.1.2. Runge-Kutta de quarta ordem

O método de Runge-Kutta de quarta ordem, ou simplesmente RK4, é um dos mais utilizados para obter soluções aproximadas de valor inicial, ele é obtido a partir do mesmo procedimento descrito anteriormente, porém com aproximação de quarta ordem, resultando em mais três interações na aproximação da equação (9), conforme os quatro passos do algoritmo apresentado a seguir. A ordem do método indica o número de pontos usado em um subintervalo para determinar o valor da inclinação, portanto, o método de Runge-Kutta de quarta ordem utiliza a inclinação em quatro pontos [¹⁵].

O método RK4 tem um erro de truncamento local proporcional a h⁵, ou seja, é três ordens de grandeza mais preciso do que o método de Euler que tem um erro da ordem h². Após n passos, o erro acumulado é proporcional a nh⁵. A seguir apresentamos as quatro interações presentes no RK4:

(10)$$k_1=F\mathrm{}(x_n,y_n),$$

(11)$$k_2=F\mathrm{}\left(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1\right)$$

(12)$$k_3=F\mathrm{}\left(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2\right)$$

(13)$$k_4=F\mathrm{}\left(x_n+h,y_n+h{k}_3\right)$$

(14) $$y(x_{n+1})\cong y(x_n)+h\mathrm{}\left(\frac{k_1+2k_2+2k_3+k_4}{6}\right).$$

Com o objetivo de comparar a eficiências dos dois métodos resolveremos a seguinte equação diferencial ordinária

(15)$$y^{\prime }(x)-y(x)=0,$$

com condição inicial y(0) = 1. É fácil verificar que a solução analítica dessa simples EDO é a função exponencial y(x) = e^x. Na Figura (1) são mostrados os gráficos dessa função a partir da condição inicial anterior pelo método analítico (em vermelho), pelo método de Euler (em verde) e pelo método de Runge-Kutta de quarta ordem (em azul) os dois últimos usando o tamanho do passo h = 0.1. Nota-se que o erro cumulativo do método de Euler mencionado acima se torna visível para uma grande quantidade de passos, enquanto que a solução do método Runge-Kutta de quarta ordem permanece com uma boa aproximação, mesmo com grandes quantidades de passos.

Figura 1
Gráfico analítico da função exponencial em vermelho para comparação com as soluções obtidas pelo método de Euler (pontos verdes) e Runge-Kutta de quarta ordem (pontos azuis).

O conjunto de equações apresentado na equação (10) descreve a implementação do método de Runge-Kutta de quarta ordem para uma única equação diferencial, ou ainda, uma equação diferencial de primeira ordem. O método pode ser utilizado para resolver equações diferenciais de segunda ordem se transformarmos em um sistema com duas equações diferenciais de primeira ordem a partir da introdução da variável adicional z = dy/dx.

Para o pêndulo amortecido, por exemplo, observamos equações diferenciais de segunda ordem. O mesmo procedimento pode ser utilizado para analisar outros problemas que envolvem sistemas de equações diferenciais acopladas, como o estudo de modelos epidemiológicos [¹⁶,¹⁷,¹⁸], por exemplo. Agora, faremos a descrição matemática do pêndulo amortecido, nosso sistema de estudo.

Os códigos da implementação foram escritos em Python e serão apresentados nos apêndices. Nas seções seguintes, montamos os diagramas de fase do pêndulo amortecido com diferentes tipos de amortecimentos usando o método escrito e também deixamos disponíveis seus códigos nos apêndices.

3. Pêndulo com Amortecimento

Estamos interessados em descrever o pêndulo simples constituído por um fio inextensível de comprimento l com uma massa m em sua extremidade movendo-se num fluido viscoso com uma força resistiva $\overrightarrow{f_r}$ proporcional à velocidade $\overrightarrow{\nu }$,¹

(16)$$\overrightarrow{f_r}=-b\overrightarrow{\nu }=-b(l\dot{\theta })\widehat{\theta }.$$

com b a constante de amortecimento e o sinal menos indicando que a força é contrária ao vetor velocidade. Encontramos a equação de movimento do pêndulo a partir da segunda lei de Newton, que será dada pela seguinte equação diferencial de segunda ordem

(17)$$m(l\ddot{\theta })=-b(l\dot{\theta })-mg\text{sen}\theta \Rightarrow \ddot{\theta }+\frac{b}{m}\dot{\theta }+\frac{g}{l}\text{sen}\theta =0,$$

ou ainda

(18)$$\ddot{\theta }+\gamma \dot{\theta }+{\omega }_0^2\mathrm{}\mathrm{sen}\theta =0.$$

onde γ = b/m e ${\omega }_0=\sqrt{g/l}$. Na aproximação de pequenas oscilações, a equação (18) é escrita na forma

(19)$$\ddot{\theta }+\gamma \dot{\theta }+{\omega }_0^2\theta =0.$$

Como era de se esperar a aplicação da lei de Newton resultou em uma equação diferencial de segunda ordem. Para uma dada equação diferencial, cada par de condição inicial descreve uma dinâmica única do sistema estudado. A evolução temporal de um sistema pode ser descrita a partir de suas condições iniciais (posição e velocidade).

Dependendo da relação entre γ e ω₀ podemos ter três tipos de soluções:

1. $\gamma <2{\omega }_0$: amortecimento subcrítico;

2. $\gamma =2{\omega }_0$: amortecimento crítico;

3. $\gamma >2{\omega }_0$: amortecimento supercrítico.

Com isso, podemos estudar os diferentes tipos de amortecimentos analisando γ. Para nos concentrar apenas no resultado da variação de γ em cada caso, mostraremos os resultados dos três amortecimentos com as mesmas condições iniciais, a amplitude angular inicial θ₀ = 0,01 rad e a velocidade angular inicial ${\dot{\theta }}_0=0,0\text{rad/s}$. Note que, para θ₀ = 0,01, a diferença entre θ₀ e sen(θ₀) é de menos de 0,002%. Em todos os casos, escolhemos g = 9,81m/s² e l = 1,5 m, consequentemente ω₀ é o mesmo. O estado do sistema é completamente determinado a partir de duas variáveis, a posição θ(t) e a velocidade $\dot{\theta }(t)$.

Pode-se, então, estudar a dinâmica em um espaço formado por essas variáveis, o espaço de fase, um espaço 2n dimensional para um sistema com n graus de liberdade. Assim, no sistema que estamos estudando, se resume a um plano. No espaço de fase, o estado varia ao longo do tempo, descrevendo uma trajetória. Cada par de condições iniciais descreve uma única trajetória e o conjunto de todas as trajetórias possíveis para todas as condições iniciais forma o diagrama de fase. Assim, a análise de trajetórias do diagrama de fase pode auxiliar o estudo físico de determinado sistema. Por exemplo, em sistemas conservativos, as trajetórias serão curvas fechadas enquanto em sistemas dissipativos, as curvas irão espiralar e tenderão à origem com o passar do tempo. Esse método se torna especialmente interessante no estudo da estabilidade de pontos de equilíbrio e na análise de sistemas não-lineares. No sistema que estamos lidando, esperamos observar esse comportamento, ou seja, trajetórias que terminam no ponto $(\theta ,\dot{\theta })=(0,0)$.

Para facilitar a visualização, os gráficos apresentados nas seções seguintes apresentam resultados limitados a t = 15s, pois os pontos gerados pela solução numérica se tornam bastante próximos dos resultados analíticos, praticamente indistinguíveis a olho nu.

3.1. Amortecimento subcrítico

A solução analítica para a equação diferencial (19) no caso em que γ < 2ω₀ é dada por [¹⁹]

(20)$$\theta (t)=A{e}^{-\frac{\gamma }{2}t}\cos (\omega t+\varphi ),$$

com $\omega =\sqrt{{\omega }_0^2-\frac{{\gamma }^2}{4}}$ e as constantes A e ϕ são determinados pelas condições iniciais, como dissemos anteriormente, θ₀ = θ(0) = 0,01 rad e ${\dot{\theta }}_0=\dot{\theta }(0)=0,0\text{rad/s}$,

(21)$$\begin{array}{lll}A\hfill & =\hfill & \frac{{\theta }_0}{\cos \left(\varphi \right)},\hfill \\ \varphi \hfill & =\hfill & \arctan \mathrm{}\left[-\frac{1}{\omega }\left(\frac{{\dot{\theta }}_0}{{\theta }_0}+\frac{\gamma }{2}\right)\right],\hfill \end{array}$$

onde escrevemos A em função de ϕ para simplificar a equação. A velocidade angular é obtida derivando a equação (20) no tempo, com isso obtemos

(22)$$\dot{\theta }(t)=A{e}^{-\frac{\gamma }{2}t}\left[\omega \text{sen(}\omega t+\varphi \text{)+}\frac{\gamma }{2}\cos (\omega t+\varphi )\right].$$

Na Figura 2 a seguir são apresentados os gráficos das soluções analítica e numérica para a posição angular e a velocidade angular usando como condições iniciais. No apêndice A apresentamos o código em Python do Runge-Kutta para este caso do pêndulo subcrítico. Observamos na Figura (2) que a amplitude das oscilações diminui com o tempo, característica do envelope exponencial que multiplica a função senoidal na equação (20). Devido ao amortecimento, a amplitude irá diminuir até que o pêndulo pare devido à perda de energia.

Figura 2
Nas figuras (a) e (b) vemos os gráficos das funções θ₀ rad e $\dot{\theta }(t)$ descritas no texto, obtidas para o pêndulo com amortecimento subcrítico. Os gráficos foram obtidos para γ = 0,1ω_o. Os pontos azuis representam as soluções numéricas usando o método de RK enquanto a linha contínua vermelha representa as soluções analíticas.

A trajetória resultante no espaço de fases, construída a partir das soluções para θ(t) e $\dot{\theta }(t)$ é mostrado na Figura (3). Como descrito anteriormente, a trajetória é uma curva em espiral partindo do ponto dado pelas condições iniciais e seguindo para a origem.

Figura 3
Trajetória no espaço de fases obtida para o pêndulo com amortecimento subcrítico com condições iniciais θ₀ = 0,01 rad e ${\dot{\theta }}_0=0,0\text{rad/s}$. A linha contínua vermelha corresponde ao resultado analítico enquanto os pontos azuis foram obtidos pela solução numérica.

3.2. Amortecimento crítico

Se, na equação (19), admitirmos que a constante $b=2m\sqrt{g/l}$ temos o caso do pêndulo com amortecimento crítico. Neste caso, a constante γ é dada por γ = 2ω₀ e a solução da equação diferencial de segunda ordem é dada pela equação [¹⁹]

(23)$$\theta (t)=(A+Bt)e^{-\frac{\gamma }{2}t}.$$

a constante A é obtida a partir da condição inicial θ(0) = θ₀, o que resulta em A = θ₀. A constante B é determinada pela segunda condição inicial $\dot{\theta }(0)={\dot{\theta }}_0$, o que fornece $B={\dot{\theta }}_0+{\omega }_0{\theta }_0$. Com essas constantes determinadas reescrevemos a equação (23) como:

(24)$$\theta (t)=[{\theta }_0+({\dot{\theta }}_0+{\omega }_0{\theta }_0)t]e^{-\frac{\gamma }{2}t},$$

e após a derivação a equação para a velocidade angular fica,

(25)$$\dot{\theta }(t)=[{\dot{\theta }}_0-{\omega }_0({\dot{\theta }}_0+{\omega }_0{\theta }_0)t]e^{-\frac{\gamma }{2}t}.$$

Na Figura (4) mostramos o comportamento das funções θ(t) da equação (24) e $\dot{\theta }(t)$ da equação (25). A equação (24) para o pêndulo com amortecido crítico mostra que o valor da variável angular diminui rapidamente com o tempo devido ao fator exponencial. Esse comportamento pode ser observado na Figura 4(a), construída utilizando o método de Runge-Kutta, cujo código está no apêndice B. Por outro lado, a velocidade angular pode apresentar uma mudança de comportamento devido ao sinal negativo observado na equação (25) e apresentado na Figura 4(b). A velocidade angular diminui rapidamente o seu valor até atingir um mínimo, quando volta a crescer até se anular.

Figura 4
Na figura, vemos o gráfico das funções (a) θ(t) e (b) $\dot{\theta }(t)$ para o pêndulo com amortecimento crítico, descrito no texto. Os gráficos foram obtidos para γ = 2ω_o. Em ambos os gráficos, um ponto preto maior que os demais assinala as condições iniciais. Os pontos azuis representam as soluções numéricas usando o método de RK enquanto a linha contínua vermelha representa as soluções analíticas.

Com o método de Runge-Kutta obtemos os gráficos apresentados na Figura (4), que mostram como a posição angular e a velocidade angular variam com o tempo para o pêndulo com amortecimento crítico, e a trajetória no diagrama de fase obtida para esse pêndulo é mostrada na Figura (5). Percebemos que inicialmente, devido ao fator exponencial, a partícula acelera rapidamente e depois desacelera para tempos maiores uma vez que a exponencial decresce mais rapidamente até atingir a velocidade angular nula. Diferente do resultado anterior, não vemos uma espiral. A trajetória apresenta uma concavidade, pois a velocidade aumenta em módulo, mas com sinal negativo. Então, a velocidade passa a diminuir seu módulo até chegar a zero.

Figura 5
Trajetória no espaço de fases obtida para o pêndulo com amortecimento crítico com condições iniciais θ₀ = 0,01 rad e ${\dot{\theta }}_0=0,0\text{rad/s}$ (apresentadas no gráfico como um ponto preto maior que os demais). Os pontos azuis foram obtidos a partir da solução numérica usando o método de RK enquanto a linha contínua vermelha representa a solução analítica.

3.3. Amortecimento supercrítico

O amortecimento supercrítico do pêndulo é observado quando a constante γ > 2ω₀. A solução da equação (19) é [¹⁹]

(26)$$\theta (t)=(A{e}^{\omega t}+B{e}^{-\omega t})e^{-\frac{\gamma }{2}t},$$

onde

(27)$$\omega =\sqrt{{\left(\frac{\gamma }{2}\right)}^2-{\omega }_0^2},$$

(28)$$A=\frac{1}{2}\left[{\theta }_0+\frac{1}{\omega }\left({\dot{\theta }}_0+\frac{\gamma {\theta }_{\text{0}}}{2}\right)\right],$$

(29)$$B=\frac{1}{2}\left[{\theta }_0-\frac{1}{\omega }\left({\dot{\theta }}_0+\frac{\gamma {\theta }_{\text{0}}}{2}\right)\right],$$

A velocidade é obtida derivando a equação equação (26), portanto

(30)$$\dot{\theta }(t)=\left[A\left(\omega -\frac{\gamma }{2}\right)\mathrm{}{e}^{2\omega t}-B\left(\omega +\frac{\gamma }{2}\right)\right]\mathrm{}{e}^{-\left(\omega +\frac{\gamma }{2}\right)t}.$$

Para verificar o comportamento das expressões anteriores utilizaremos o valor γ = 5ω₀. Os gráficos da posição angular e velocidade angular para o pêndulo com amortecimento supercrítico estão mostrados na Figura (6) enquanto o código do Runge-Kutta está no apêndice C. Como o caso anterior, a variável angular vai rapidamente a zero enquanto a velocidade angular começa assumindo valores negativos antes de diminuir seu módulo até chegar a zero. Contudo, para o pêndulo com amortecimento supercrítico, observamos que o gráfico apresenta uma curva menos aguda do que o caso anterior, enquanto a velocidade atinge um mínimo menor, em módulo, do que o anterior.

Figura 6
Nas figuras (a) e (b), são mostradas o gráfico das funções θ(t) e $\dot{\theta }(t)$ para o pêndulo com amortecimento supercrítico com γ = 5ω₀, respectivamente. A solução analítica é apresentada como uma linha contínua vermelha e as soluções obtidas com o método de RK como pontos azuis nas imagens. Em ambos os gráficos, um ponto preto maior que os demais assinala as condições iniciais θ₀ = 0,01 rad e ${\dot{\theta }}_0=0,0\text{rad/s}$.

A trajetória no espaço de fase para o pêndulo com amortecimento supercrítico é mostrada na Figura (7). Nesse caso, a trajetória tem uma aparência mais aguda que o anterior, em virtude dos menores valores atingidos pela velocidade. Novamente, não vemos espirais nessa imagem devido ao amortecimento mais intenso. A trajetória rapidamente chega à origem.

Figura 7
A trajetória no espaço de fase obtida com condições iniciais θ₀ = 0,01 rad e ${\dot{\theta }}_0=0,0\text{rad/s}$ (mostrada no gráfico como um grande ponto preto) é apresentada no gráfico em vermelho (solução analítica) e como pontos azuis (solução numérica).

4. Regimes de Validade

Nas seções anteriores, mostramos as soluções para o pêndulo amortecido em regimes cujas equações diferenciais apresentam solução analítica. Essas soluções são válidas quando sen(θ) ≈ θ. Para mostrar as soluções usando o método de Runge-Kutta, escolhemos θ = 0,01 rad, o que está condizente com a aproximação. Para ilustrar, o gráfico na Figura 8a abaixo mostra a comparação entre entre valores exatos de sen(θ) e os valores de θ em radianos, para o ângulo variando entre 0.0 e 1.5 radiano com pontos distanciados de 0.02 radianos. No inset destacado no gráfico está indicado um zoom na região para ângulos menores do que 0.30 radianos, no qual observa-se uma concordância muito boa entre as duas funções.

Figura 8
Comparação entre as funções sen(θ), pontos azuis, e θ, a linha contínua vermelha, mostrada na figura (a) evidencia o regime no qual a aproximação seno(θ) ≈ θ. é válida. O inset da figura mostra as mesmas funções para ângulos menores do que 0,30 radianos. Em (b), mostramos o erro percentual entre as duas funções. O inset de (b) mostra o erro percentual para ângulos menores do que 0,30 radianos.

Na Figura 8b, mostramos o aumento do erro percentual entre os valores de sen(θ) e os valores de θ, calculado como |sen(θ) − θ|/sen(θ). Nota-se, no inset, que o erro da aproximação é de menos de 2% para ângulos menores do que 0,30 radianos.

Escolhemos mostrar nossos resultados no regime de pequenas oscilações pois as equações diferenciais possuem solução analítica, de forma que podemos comparar o resultado exato com o resultado obtido pelo método numérico. Uma vez que observamos a precisão do método, é interessante compará-lo com resultados onde a aproximação deixa de ser válida. Na Figura 9 mostramos três gráficos contendo as trajetórias no espaço de fase no regime criticamente amortecido.

Figura 9
Trajetórias no espaço de fase para o sistema criticamente amortecido no regime de pequenas oscilações (linha contínua vermelha) e usando o Runge-Kutta de quarta ordem (pontos azuis) com a posição angular inicial de 0,3 rad (a), 0,6 rad (b) e 0,9 rad (c).

A Figura 9a mostra o caso em que a posição angular inicial é 0,3 radianos. Como discutimos, o erro na aproximação do seno pelo ângulo é de menos de 2%. Esse resultado se reflete no gráfico que apresenta uma boa concordância entre os dois métodos. Contudo, aumentando θ₀ para 0,6 radianos (Figura 9b) e 0,9 radianos (Figura 9c), as soluções se distanciam perceptivelmente. Nesses casos, a solução analítica discutida anteriormente deixa de representar satisfatoriamente o sistema apresentado e é preciso usar métodos numéricos para resolver as equações diferenciais. Resultados semelhantes são mostrados na Figura 10 para os regimes subamortecido 10(a) e sobreamortecido 10(b). Observamos que, à medida que o ângulo inicial aumenta, é bastante difícil perceber a diferença entre as duas soluções no caso subamortecido. Contudo, as trajetórias para o regime sobreamortecido já revelam uma diferença considerável entre a solução analítica e a solução numérica, assim como visto no caso criticamente amortecido.

Figura 10
Comparação entre as trajetórias no espaço de fase para os regimes subamortecido (a) sobreamortecido (b). As soluções obtidas no limite de pequenas oscilações são apresentadas como linhas contínuas vermelha e os pontos azuis representam o resultado obtido usando o Runge-Kutta de quarta ordem na equação diferencial sem aproximação. A posição angular inicial aumenta de cima para baixo, sendo a primeira linha 0,3 radianos, a segunda linha 0,6 radianos e a última linha 0,9 radianos.

5. Conclusões

Diferentes métodos numéricos são usados por físicos a depender do problema estudado pelo pesquisador. Um dos mais poderosos métodos para resolver equações diferenciais de primeira ordem com problema de valor inicial é o método de Runge-Kutta de quarta ordem. Nesse trabalho, discutimos esse método e sua implementação em linguagens de computação (escolhemos Python para isso). Mostramos como implementar o método para resolver um problema tradicional da física, o pêndulo amortecido. Confirmamos os resultados ao compará-los com as soluções analíticas desses sistemas físicos. Esperamos que nosso trabalho sirva de ajuda como um estudo introdutório para estudantes que realizarão projetos de pequisa que necessitem obter soluções numéricas de equações diferenciais.

Agradecimentos

Os autores são extremamente gratos ao Instituto Federal do Sertão Pernambucano e ao árbitro cujos comentários foram essenciais para aprimorar e enriquecer o texto apresentado nesse trabalho.

Material Suplementar

O seguinte material suplementar está disponível online

Apêndice A Código: Pêndulo com amortecimento subcrítico.

Apêndice B Código: Pêndulo com amortecimento crítico.

Apêndice C Código: Pêndulo com amortecimento supercrítico.

Referências

• ^[1]
  W. Kinzel e G. Reents, Physics by Computer: Programming Physical Problems Using Mathematica® and C (Springer, Berlin, 1997).
• ^[2]
  M.R. Matthews, Science & Education 13, 689 (2004).
• ^[3]
  L.P. Fulcher, B.F. Davis, American Journal of Physics 44, 51 (1976).
• ^[4]
  L.N. Gonçalves, FernandesJ., FerrazA., SilvaA.G. e P.J. Sebastião, American Journal of Physics 88, 962 (2020).
• ^[5]
  R.D. Peters, ArXiv preprint physics/0306081 (2003).
• ^[6]
  M. Kavyanpoor e S. Shokrollahi, Journal of King Saud University-Science 31, 14 (2019).
• ^[7]
  D.X. Andrade, P.H.R. Anjos e P.E.G. Assis, Revista Brasileira de Ensino de Física 39, e1307 (2017).
• ^[8]
  F.M.S. Lima, European journal of physics 29, 1091 (2008).
• ^[9]
  R.N. Suave e J.A Nogueira, Revista Brasileira de Ensino de Física 38, 2501 (2016).
• ^[10]
  J.C. Fernandes, P.J. Sebastião, L.N. Gonçalves e A. Ferraz, European Journal of Physics 38, 045004 (2017).
• ^[11]
  K. Johannessen, European Journal of Physics 35, 035014 (2014).
• ^[12]
  T.A. Sobral, De V.H. Holanda, F.C.B. Leal e T.T. Saraiva, Journal of Physics D: Applied Physics 54, 255308 (2021).
• ^[13]
  M. Pigneur e J. Schmiedmayer, Physical Review A 98, 063632 (2018).
• ^[14]
  G.B. Arfken e H.J. Weber, Física Matemática: Métodos matemáticos para engenharia e física (Elsevier, 2007).
• ^[15]
  A. Gilat e V. Subramanaim, Métodos numéricos para engenheiros e cientistas: uma introdução com aplicações usando o MATLAB (Bookman, Porto Alegre, 2008).
• ^[16]
  M. Al-Raeei, Modeling Earth Systems and Environment 8, 1443 (2022).
• ^[17]
  A.J. Arenas, J.A Moraão e J.C. Cortês, Computers Mathematics with Applications 56, 670 (2008).
• ^[18]
  O.O. Onyejekwe, A. Tigabie, B. Ambachew e A. Alemu, Journal of Applied Mathematics and Physics 7, 148, (2019).
• ^[19]
  H.M. Nussenzveig, Curso de Física Básica (Edgard Blucher, Sãoo Paulo, 2002), v. 2, 4 ed.

• 1
  Em coordenadas polares: $\overrightarrow{v}=\dot{r}\widehat{r}+r\dot{\theta }\widehat{\theta }$. Quando r = l = const., teremos $\overrightarrow{v}=(l\dot{\theta })\widehat{\theta }$.

Datas de Publicação

• Publicação nesta coleção
  10 Jun 2024
• Data do Fascículo
  2024

Histórico

• Recebido
  12 Mar 2024
• Revisado
  15 Abr 2024
• Aceito
  17 Abr 2024