Análise numérica e visualização de osciladores via Python

Numerical analysis and visualization of oscillators via Python

Resumo

Apresentamos nesse trabalho o uso da linguagem de programação Python e suas bibliotecas VPython, NumPy, Pandas e Matplotlib como recursos de desenvolvimento de simulações computacionais relacionadas a problemas de osciladores para serem aplicados nas aulas de física nos cursos de graduação em física, exatas e engenharias como uma alternativa de experimento virtual. A discussão foi feita no comparativo entre os resultados provenientes de soluções exatas e soluções aproximadas obtidas através de métodos numéricos e na construção de simuladores visuais em três tipos de osciladores: sistema massa-mola com rotação, pêndulo físico acoplado a uma mola e um sistema de oscilador acoplado. O estudo comparativo, entre as soluções exatas e as soluções aproximadas, foi feito usando os métodos numéricos de Euler, de Heun e o método de Runge-Kutta de quarta ordem (RK4), onde foi considerado três passos de tempo, $\delta t=t_{n+1}-t_n$, de 0.2 e 0.05 para cada um dos métodos. Os resultados mostraram que o método de RK4 forneceu uma melhor aproximação da solução exata nos valores de posição e velocidade dos osciladores. A construção e visualização da animação dos osciladores foram feitas através da implementação computacional do método de Euler e o uso da biblioteca VPython, proporcionando um ambiente de aprendizado mais interativo.

Palavras-chave:
Programação Python; Métodos Numéricos; Osciladores; Simuladores Visuais

Abstract

In this work, we present the use of the Python programming language and its libraries VPython, NumPy, Pandas and Matplotlib as resources for developing computational simulations related to oscillator problems to be applied in physics classes in undergraduate courses in physics, exact sciences and engineering as an alternative to virtual experiments. The discussion is made by comparing the results from exact solutions and approximate solutions obtained through numerical methods and by building visual simulators in three types of oscillators: mass-spring system with rotation, physical pendulum coupled to a spring and a coupled oscillator system. The comparative study, between the exact solutions and the approximated solutions, was done using the numerical methods of Euler, Heun and the fourth-order Runge-Kutta method (RK4), where three time steps were considered, $\delta t=t_{n+1}-t_n$, of 0.2 and 0.05. The results showed that the RK4 method presented a better approximation of the exact solution in the position and velocity values of the oscillators. The construction and visualization of the oscillator animation were done through the computational implementation of the Euler method and the use of the VPython library, providing a more interactive learning environment.

Keywords
Python Programming; Numerical Methods; Oscillators; Visual Simulators

1. Introdução

Do ponto de vista experimental, a visualização de eventos reais realizados em laboratório de física é algo que nem sempre é fácil de perceber, uma vez que tais experimentos ocorrem de maneira muito rápida, em intervalos de tempos tão curtos que, algumas vezes, estão aquém de serem observáveis pela percepção humana (da audição, visão, etc).

Relacionar os conceitos de física com o cotidiano nem sempre é uma tarefa fácil. A dificuldade de aprendizagem desta ciência é mostrada nos altos índices de desistência e reprovação nos vários níveis de ensino [¹, ², ³, ⁴]. No entanto, é possível aproximar o aluno do conhecimento através de atividades experimentais virtuais que favoreçam o desenvolvimento do espírito crítico, da elaboração de hipóteses, raciocínio lógico e uma melhor percepção da relação ciência e tecnologia [⁵].

Visando proporcionar um entendimento mais completo do fenômeno físico a partir de atividades experimentais, a utilização de novas tecnologias, tais como simulações computacionais, como ferramenta de clarificação de fenômenos físicos podem auxiliar com mais rapidez o entendimento do fenômeno de interesse. Simulações e animações oferecem ilimitado potencial de aplicação [⁶], o que lhes permite serem chamados de laboratórios virtuais, podendo incentivar e motivar a busca pelo entendimento de princípios teóricos das ciências naturais. A utilização das tecnologias computacionais também possibilita a criação de ambientes virtuais de aprendizagem, nos quais os estudantes podem interagir e colaborar de maneira ativa na exploração dos conceitos físicos.

A ferramenta computacional tem um papel muito importante para a pesquisa na área de física. Existe uma série de fenômenos físicos que são difíceis de serem reproduzidos em um laboratório com experimentos reais e também são extremamente complicados de serem solucionados do ponto de vista da matemática. Diante dessas dificuldades, a utilização de alguma liguagem de programação na computação científica se apresenta como um importante instrumento auxiliar no estudo de fenômenos físicos, pois ela serve como um recurso valioso para investigar e compreender esses fenômenos de maneira mais eficaz. A adoção dessa abordagem, além de facilitar a modelagem e simulação, também amplia as possibilidades de análise e interpretação dos resultados.

Vários pesquisadores sugerem que atividades mediadas por recursos tecnológicos e multimídias sejam aplicadas em sala de aula [⁷, ⁸, ⁹, ¹⁰, ¹¹]. A utilização de uma ferramenta computacional nos permite ter uma análise mais completa de um determinado fenômeno. No estudo de oscilação não amortecida e amortecida, por exemplo, é possível obter de maneira rápida informações acerca do seu movimento através de gráficos do deslocamento, velocidade e energia por meio de análise numéricas, como feito nas Refs.[¹², ¹³, ¹⁴, ¹⁵, ¹⁶, ¹⁷, ¹⁸].

Uma linguagem de programação que é usada como referência no campo da física é o Python. Esta linguagem conta com pacotes para análise e manipulação de dados com a biblioteca Pandas, além de pacotes avançados matemáticos, o NumPy, incluindo funcionalidades gráficas como a Matplotlib, que auxiliam na criação de códigos para, por exemplo, resolver numericamente problemas de equações diferenciais [¹⁹]. O Python é apresentado como uma linguagem de programação fácil de aprender e robusta, o que se torna apropriada para aqueles com pouca ou nenhuma experiência em programação além de ser muito adequada para criar modelos interativos de sistemas físicos em 3D utilizando a biblioteca VPython [²⁰]. Outras biliotecas como o Sympy e o Manim também são utlizadas respectivamente para manipulação simbólica, a qual se apresenta como uma ferramenta alternativa a sistemas como Mathematica ou Maple, e construções de animações gráficas, projetada para criar animações de alta qualidade para matemática, ciências e tópicos relacionados a dados. Apresentamos na Tab. 1 um resumo dessas bibliotecas e suas especificidades.

Tabela 1
Bibliotecas do Python.

Biblioteca	Especificidade
Pandas	análise e manipulação de dados
Numpy	manipulação de computação numérica
Matiplotlib	construções de gráficos
VPython	construções de modelos interativos 3D
Sympy	manipulação simbólica
Manim	construções de animações gráficas 2D

Este trabalho tem como finalidade apresentar o uso da linguagem de programação Python e suas bibliotecas VPython, NumPy, Pandas e Matplotlib como recursos de desenvolvimento de simulações computacionais relacionadas a problemas de osciladores e para serem aplicados nas aulas de física nos cursos de graduação em física, exatas e engenharias como uma alternativa de experimento virtual. A discussão se baseia na construção de simuladores visuais e no comparativo entre os resultados provenientes de soluções exatas de equações diferenciais de segunda ordem e soluções aproximadas obtidas através de métodos numéricos em três tipos de osciladores: sistema massa-mola com rotação, pêndulo físico acoplado a uma mola e um sistema de oscilador acoplado. Para a implementação computacional desses osciladores, os métodos numéricos utilizados foram o método de Euler, o método de Heun e o método de Runge-Kutta de quarta ordem.

O método de Euler foi introduzido por Leonhard Paul Euler em seu livro em 1768 [²¹]. Esse método, também conhecido como método da reta tangente, é utilizado na resolução aproximada de equações diferenciais ordinárias de primeira ordem com um problema de valor incial dado, utilizando a descrição da série de Taylor [²², ¹⁶].

O método de Heun, também conhecido como método Preditor-Corretor, é uma versão aperfeiçoada do método de Euler [²³]. O método, apresenta uma abordagem para melhorar a estimativa da inclinação por meio da utilização de duas derivadas, uma no ponto inicial e outra no ponto final, resultando em uma melhor aproximação em comparação ao valor real.

Os métodos de Runge-Kutta foram inicialmente desenvolvidos por C. Runge em 1895 [²⁴], enquanto W. Kutta estendeu esta classe de métodos em 1901 [²⁵]. Os métodos representam toda uma família de algoritmos para a solução numérica de problemas de valor inicial para EDOs de primeira ordem ou sistemas destes fundamentados na comparação de um polinômio de Taylor apropriado para eliminar o cálculo das derivadas, efetuando diversas avaliações da função a cada passo [²⁶, ²⁷]. Dentre esses métodos, o de quarta ordem (RK4) é um dos mais usados para obter soluções numéricas aproximadas. Embora seja mais complexo do que o método de Euler e o método de Heun, o método RK4 é também mais preciso em seus resultados.

O ponto chave desenvolvido nesse trabalho integra Física, métodos numéricos e física computacional, o qual relaciona diretamente o conhecimento tradicional às novas técnicas de análise e estudo em física. Diante disso, o artigo foi organizado da seguinte forma: na próxima seção, são apresentados os modelos de osciladores utilizados para o estudo. Na seção 3, apresentamos os métodos numéricos utilizados na investigação. Na seção 4 os resultados e discussões acerca dos métodos utilizados e a construção dos simuladores visuais para os diferentes modelos de osciladores e por fim, na seção 5, as conclusões obtidas.

2. Osciladores

Quando o movimento se repete de maneira regular no decorrer do tempo, é considerado periódico, e o intervalo entre duas ocorrências equivalentes é denominado período do movimento. Os sistemas físicos que realizam esse tipo de movimento são chamados deosciladores.

Para uma partícula que executa um movimento periódico, a força resultante que atua sobre ela é nula na posição de equilíbrio. Durante a oscilação, a posição da partícula, a sua velocidade e a aceleração variam periodicamente com o tempo e, consequentemente, a força que atua sobre ela também varia. A descrição matemática desse tipo de movimento envolve o uso de funções trigonométricas senos e cossenos ou de combinações dessas funções, conhecidas como funções harmônicas, o que leva o movimento periódico a ser conhecido como movimento harmônico.

O movimento oscilatório é um importante tópico de estudo na física, devido à frequencia que este tipo de evento pode ocorrer, desde situações simples como o movimento de uma cadeira de balanço, como em situações mais sofisticadas, como em circuitos elétricos oscilantes RLC [²⁸, ²⁹] e na emissão de sons em instrumentos musicais [³⁰, ³¹], por exemplo. Os mais variados estudos de sistemas oscilatórios foram feitos usando simulações computacionais e cálculo numérico apartir da linguagem de programação Python e sua biblioteca gráfica VPython. Entre eles, temos a análise do movimento do pêndulo simples sob a influência de uma força externa [³²]; a visão geral do uso do VPython para modelar e simular alguns fenômenos simples na física com visualização explícita dos comportamentos desses fenômenos sob condições variadas [³³]; o estudo de um pêndulo harmônico como exemplo didático de integração de laboratório construído com peças LEGO e implemetação númérica Python [³⁴]; abordagem computacional para resolver o problema do oscilador harmônico simples usando o método de Runge-Kutta [³⁵]; investigação da relação entre sincronização de pêndulo simples e duplo e sincronização do mundo real implementado em Python [³⁶]; e o estudo do caos determinístico através da simulação do movimento de um pêndulo duplo, permitindo a visualização de trajetórias intrincadas e a análise de padrões emergentes com as linguagens de programação Python e Fortran [³⁷].

Nesta seção, apresentamos e discutimos os três modelos de osciladores utilizados no trabalho: sistema massa-mola com rotação, pêndulo físico acoplado a uma mola e um oscilador harmônico acoplado. Trataremos aqui de suas principais características e análise acerca de suas equações de movimento.

2.1. Modelo 1: sistema massa-mola com rotação

Um oscilador harmônico pode ser representado por mola de constante elástica k, de massa desprezível, que tem uma extremidade fixada numa parede e a outra, presa a um disco de massa m e raio R, que rola sem deslizar numa superfície horizontal, como mostrado na Fig. 1. Nesta situação, existem duas forças relevantes que atuam sobre o corpo: uma é resultado do estiramento ou compressão da mola e a outra é devido ao atrito na superfície, causando o movimento de rolamento do disco. Quando a mola está relaxada, a força resultante sobre o disco é nula e o sistema está em equilíbrio.

Figura 1
Oscilação de um disco preso a uma mola sem deslizamento.

O disco se move em movimento harmônico simples quando é deslocado de sua posição de equilíbrio e depois liberado, seguindo uma trajetória retilínea unidirecional e executando um movimento periódico de translação e rotação. A força que atua na mola é uma força restauradora linear, lei de Hooke, dada pelo produto entre a constante elástica k da mola e o deslocamento x do bloco em torno de sua posição de equilíbrio, $\overrightarrow{F}=-k\overrightarrow{x}$. Sabemos também que a Segunda Lei do Movimento de Newton nos diz que $\sum \overrightarrow{F}=m\overrightarrow{a}$, sendo assim

(1)$$ma=-kx-f_{atrito}.$$

Note que, para o caso da força de atrito, f_atrito, ser nula teremos, na equação Eq01, o problema físico bem conhecido que é o oscilador harmônico simples. Portanto, nesse caso, o disco executará apenas um movimento periódico de translação, pois a força de atrito é responsável pelo rolamento do disco. A força de atrito estatico é obtida aplicando a segunda na lei na forma angular, $\sum \overrightarrow{\tau }=I\overrightarrow{\alpha }$. Sendo assim, temos que o torque, $\overrightarrow{\tau }$, aplicado no centro de massa do disco é dado por τ = Iα = Rf_atrito, onde α = a/R é a aceleração angular e I = mR²/2 é o momento de inércia do disco, então f_atrito = ma/2. Portanto, a equação Eq01 assume a seguinte forma

(2)$$a(t)=-\frac{2k}{3m}x(t)\Rightarrow \frac{d^{2}x(t)}{d{t}^2}=-\frac{2k}{3m}x(t),$$

que é a equação diferencial para o sistema cuja solução pode ser escrita na seguinte forma

(3)$$x(t)=x_m\text{cos}({\omega }_{0}t).$$

Essa é a equação de posição x(t) para corpo em relação a posição de equilíbrio, na qual x_m é a amplitude do movimento e ω₀ é a frequência angular,

(4)$${\omega }_0=\sqrt{\frac{2k}{3m}}.$$

A velocidade v(t) é obtida por v(t) = dx(t)/dt, assim

(5)$$v(t)=-x_m{\omega }_0\text{sen}({\omega }_{0}t).$$

2.2. Modelo 2: pêndulo físico acopladoa uma mola

Um pêndulo físico é caracterizado por oscilações semelhantes às de um pêndulo simples, mas não pode ser representado como uma massa pontual suspensa por um fio. Em vez disso, a distribuição de massa do objeto deve ser considerada na equação de movimento, o que é feito por meio do momento de inércia. Este parâmetro é fundamental, pois descreve a forma como a massa está distribuída em relação ao ponto de oscilação.

Considere que uma barra rígida de comprimento L e massa m_b está suspensa, podendo girar em torno do ponto O, por uma das suas extremidades. Na outra extremidade, a barra está ligada a um disco de raio R e de massa m_d. Uma mola de constante k, que está na posição relaxada quando a barra se encontra na posição vertical, é conectada à barra a uma distância d abaixo de seu ponto de suspensão. A barra é deslocada para a esquerda, até um ângulo θ com a direção vertical, e abandonada a partir do repouso, a qual descreve um arco de circunferência em um movimento harmônico simples, como mostra a Fig. 2(a).

Figura 2
Diagrama ilustrativo para um pêndulo físico acoplado a uma mola.

A barra está sujeita à força peso ${\overrightarrow{P}}_b$, aplicada no centro de massa da barra, à distância L/2 do ponto O, e também à força elástica exercida pela mola ${\overrightarrow{F}}_{el}$ aplicada a uma distância d do ponto O, na direção horizontal. Temos também a força peso ${\overrightarrow{P}}_d$, aplicada no centro de massa do disco, à distância R + L do ponto O, como mostrado na Fig. 2(b). Somente as componentes dessas forças, na direção perpendicular à barra, exercem torque. Nesta situação, temos a componente perpendicular da força peso da barra dada por m_bgsen(θ) e a componente perpendicular da força peso do disco dada por m_dgsen(θ), como mostra a Fig. 2(c). O torque resultante, τ_R, no pêndulo é dado pela contribuição do torque exercido pela componente perpendicular da força peso da barra τ_b = (L/2)m_bgsen(θ), contribuição do torque exercido pela força componente perpendicular da peso do disco τ_d = (R + L)m_dgsen(θ) e pela contribuição da força elástica aplicada pela mola τ_F = F_eladcos(θ). Então,

(6)$${\tau }_R=m_{b}g\frac{L}{2}\text{sen}(\theta )+m_{d}g(R+L)\text{sen}(\theta )+F_{ela}d\text{cos}(\theta ),$$

a força exercida pela mola, é F_el = kx, onde x é a deformação da mola e tem a forma x = hsen(θ). Assim:

(7)$${\tau }_R=\left[m_{b}g\frac{L}{2}+m_{d}g(R+L)+k{d}^2\text{cos}(\theta )\right]\text{sen}(\theta ).$$

Como o torque é restaurador, a equação do movimento será dada por τ_R = −Iα. Logo

(8)$$\alpha =-\frac{1}{I}\left[m_{b}g\frac{L}{2}+m_{d}g(R+L)+k{d}^2\text{cos}(\theta )\right]\text{sen}(\theta )$$

Para pequenas amplitudes, onde θ ˂˂ 1, temos que cos(θ) ≈ 1 e sen(θ) ≈ θ. Portanto, a equação(8) pode ser escrita como

(9)$$\alpha (t)=\frac{d^2\theta (t)}{d{t}^2}=-\left[\text{}\frac{(m_b+2m_d)gl+2m_{d}gR+2k{d}^2}{2I}\text{}\right]\text{}\theta (t)$$

que é a equação diferencial para o sistema, onde I é o momento de inércia total em relação ao ponto O dado por

(10)$$I=\frac{1}{3}{m}_b{L}^2+\frac{1}{2}{m}_d{R}^2+m_d{\left(R+L\right)}^2.$$

A equação(9) pode ser resolvida por

(11)$$\theta (t)={\theta }_M\text{cos}({\omega }_{0}t),$$

que é a equação da posição θ(t) em relação à posição de equilíbrio, ponto O, onde θ_M é a amplitude do movimento e ω₀ é a frequência natural do movimento chamada de frequência angular,

(12)$${\omega }_0=\sqrt{\frac{(m_b+2m_d)gl+2m_{d}gR+2k{d}^2}{2I}}.$$

A velocidade angular ω(t) é obtida por ω(t) = dθ(t)/dt, assim

(13)$$\omega (t)=-{\theta }_M{\omega }_0\text{sen}({\omega }_{0}t),$$

essa é a velocidade em função do tempo no qual o pêndulo físico gira em torno da posição de equilíbrio, ponto O, como mostrado na Fig. 2.

2.3. Modelo 3: oscilador acoplado

Osciladores acoplados são osciladores conectados de tal forma que a energia pode ser transferida entre eles. Considere dois corpos de massas iguais a m que podem se mover sem atrito sobre um trilho e estão presos por meio de duas molas de constantes elásticas k₁ e k₂, conforme mostrado na Fig. 3. Suponha que a posição de equilíbrio do sistema esteja em cada mola não esticada, sendo x₁ o deslocamento da primeira massa em relação ao seu equilíbrio e x₂ o deslocamento da segunda massa do seu equilíbrio. As forças ${\overrightarrow{F}}_1$ e ${\overrightarrow{F}}_2$ são forças restauradoras que obedecem a Lei de Hooke, como visto no Modelo 1. A força que o corpo 1 sente devido as ações das molas é −k₁x₁ + k₂(x₂ − x₁) e para o corpo 2 temos que a força é dada por −k₂(x₂ − x₁). Estas forças resultam em duas equações diferenciais quando a segunda lei de Newton é aplicada em cada um dos corpos, −(k₁ + k₂)x₁ + k₂x₂ e ma₂ = k₂x₁ − k₂x₂. Portanto

Figura 3
Oscilação para o sistema acoplado.

(14)$$\frac{d^2{x}_1(t)}{d{t}^2}=-\frac{k_1+k_2}{m}{x}_1(t)+\frac{k_2}{m}{x}_2(t)$$

e

(15)$$\frac{d^2{x}_2(t)}{d{t}^2}=\frac{k_2}{m}{x}_1(t)-\frac{k_2}{m}{x}_2(t).$$

Estas equações são extremamentes exaustivas de se obter as soluções, o que não é a proposta desse trabalho. Os detalhes das soluções podem ser vistas no Apêndice [apdA]A. As soluções para x₁(t) e x₂(t) são dadas respectivamente por

(16)$$x_1(t)=A_1\text{cos}({\omega }_{1}t)+A_2\text{cos}({\omega }_{2}t)$$

e

(17)$$x_2(t)=B_1\text{cos}({\omega }_{1}t)+B_2\text{cos}({\omega }_{2}t),$$

com A₁, A₂, $B_1=\left(k_1+\sqrt{k_1^2+4k_2^2}\right)A_1/2k_2$ e $B_2=\left(k_1-\sqrt{k_1^2+4k_2^2}\right)A_2/2k_2$ sendo as amplitudes do movimento. Para as frequências de oscilações ω₁ e ω₂ temos respectivamente

(18)$${\omega }_1=\pm {\left[\frac{k_1+2k_2}{2m}-\frac{\sqrt{k_1^2+4k_2^2}}{2m}\right]}^{1/2}$$

(19)$${\omega }_2=\pm {\left[\frac{k_1+2k_2}{2m}+\frac{\sqrt{k_1^2+4k_2^2}}{2m}\right]}^{1/2}$$

e as velocidades v₁(t) e v₂(t) são dadas por

(20)$$v_1(t)=-{\omega }_1{A}_1\text{sen}({\omega }_{1}t)-{\omega }_2{A}_2\text{sen}({\omega }_{2}t)$$

e

(21)$$v_2(t)=-{\omega }_1{B}_1\text{sen}({\omega }_{1}t)-{\omega }_2{B}_2\text{sen}({\omega }_{2}t),$$

Podemos ver ainda os vínculos do sistema através de dois modos de vibração, o simétrico e o anti-simétrico. Perceba que, nas equações(16) e (17), sempre que apenas um dos coeficientes for diferente de zero, o movimento envolverá oscilação harmônica simples (mas pertencente a ambos os corpos ao mesmo tempo). Sendo assim, para o caso de A₂ = 0 obtemos B₂ = 0 e ax₁(t) = x₂(t) = A₁cos(ω₁t), onde $a=\left(k_1+\sqrt{k_1^2+4k_2^2}\right)/2k_2$, temos um modo de vibração simétrico. Observamos então como se tivéssemos um único oscilador. Já para o caso de A₁ = 0 obtemos B₁ = 0 e bx₁(t) = x₂(t) = A₂cos(ω₂t), onde $b=\left(k_1-\sqrt{k_1^2+4k_2^2}\right)/2k_2$, temos um modo de vibração anti-simétrico, ou seja, ambos os corpos se movem em sentidos opostos. Como frequências mais altas correspondem a energias mais altas, o modo assimétrico tem uma energia mais alta.

3. Métodos Numéricos

Dizemos que uma equação diferencial ordinária (EDO) é de primeira ordem se apenas a primeira derivada da função desconhecida aparece na equação. Se quisermos encontrar uma solução para uma equação diferencial que cumpra um determinado valor inicial, dizemos que temos um problema de valor inicial do tipo

(22)$$\frac{df}{dx}=f(x,y),y(x=x_0)=y_0.$$

Para o caso de uma EDO de segunda ordem podemos reescrevê-la na forma de duas EDO’s de primeira ordem. Por exemplo, sabemos que em um problema de Física a aceleração é a segunda derivada temporal da posição, assim a = d²x/dt² e a velocidade é a primeira derivada temporal da posição v = dx/dt, então é possível estudar o sistema de interesse através do par de EDO’s de primeira ordem a = dv/dt e v = dx/dt com suas devidas condições iniciais.

Mesmo quando sabemos que o problema de valor inicial tem solução, não é necessariamente possível resolver o problema utilizando uma fórmula ou uma receita de cálculo. Sendo assim, temos que recorrer a outros meios, como métodos numéricos, para dizer algo sobre as soluções obtidas. Diante disso, discutiremos aqui três métodos numéricos para resolver sistemas físicos regidos por EDO de primeira e segunda ordem.

3.1. Método de Euler

Considere um problema que obedece uma equação diferencial ordinária f(x, y) e condição de contorno y(x = x₀) = y₀ como a equação ME:Eq01, onde y é a variável dependente e x é a variável independente. Podemos encontrar uma solução aproximada para essa equação, f(x, y), dividindo o intervalo [x₀, x_f ] em um grande número de pequenos intervalos de igual largura δx que é denominado tamanho do passo, então para n intervalos temos o passo δx = x_n+1 − x_n. Dessa forma, podemos inicialmente aplicar na equação diferencial f(x, y) os valores iniciais f(x₀, y₀) para descobrir uma primeira aproximação. Sendo assim, a estratégia do método de Euler consiste em calcular y₁ = y₀ + f(x₀, y₀)δx, que é, no primeiro passo, a aproximação y₁ da solução exata y(x₁) no ponto x₁ no ponto x₁ na equação(22). Em seguida, y₂ = y₁ + f(x₁, y₁)δx, que é, no segundo passo, a aproximação y₂ da solução exata y(x₂) no ponto x₂. Continuando esse processo sucessivas vezes chegaremos a uma fórmula geral

(23)$$x_{n+1}=x_n+\delta x$$

e

(24)$$y_{n+1}=y_n+f(x_n,y_n)\delta x$$

o qual nos dá uma solução aproximada da solução exata para a equação (22). Esse procedimento é chamado método de Euler. A Fig. 4 mostra um fluxograma que nos permite visualizar o processo feito nas interações.

Figura 4
Fluxograma método de Euler.

Temos que para uma análise numérica em um computador o processo deve seguir as seguintes instruções: (i) entrar com as condições iniciais x = x₀ e y(x₀) = y₀, o passo δx e a função f(x_n, y_n); (ii) enquanto x ˂ x_f as equações (23) e (24) são utilizadas para armazenar os dados; x ˂ x_f caso contrário, com os dados obtidos são plotados os gráficos ou guardados em uma lista (tabela de dados).

3.2. Método de Heun

Para uma melhor estimativa dos resultados numéricos, temos o método de Heun que é uma versão aperfeiçoada do método de Euler. O método consiste numa estratégia que envolve a utilização da média de duas derivadas, uma no ponto inicial e outra no ponto final. Sendo assim, modifica-se a fórmula do método Euler, equação (24), para

(25)$$y_{n+1}=y_n+\frac{1}{2}\left[f(x_n,y_n)+f(x_{n+1},y_{n+1})\right]\delta x.$$

Combinando com a equação (24), obtemos o método de Heun na forma

(26)$$y_{n+1}=y_n+\frac{1}{2}[f(x_n,y_n)+f(x_{n+1},y_n+f(x_n,y_n)\delta x)]\delta x.$$

Esse método apresenta um número maior de cálculos pois a equação (26) é calculada em cada interação garantindo então uma melhor estimativa para o resultado numérico quando comparando com a equação (24).

Uma análise numérica com a implementação do método de Heun é feita da seguinte forma em um computador: (i) entrar com as condições iniciais x = x₀ e y(x₀) = y₀, o passo δx e a função f(x_n, y_n); (ii) é calculado f(x_n, y_n); (iii) em seguida é feita a estimativa de y_n+1 usando o método de Euler; (iv) é calculado f(x_n+1, y_n+1); (v)enquanto x ˂ x_f as equações (25) e (26) são utilizadas para armazenar os dados; (vi) caso contrário, com os dados obtidos são plotados os gráficos ou guardados em uma lista (tabela de dados). A Fig. 5 apresenta um fluxograma que ilustra a implementação do método de Heun.

Figura 5
Fluxograma método de Heun.

3.3. Método de Runge-Kuta de quarta ordem

Um dos métodos numéricos mais utilizado para obter soluções aproximadas de equações diferenciais do tipo (22) são os métodos de Runge-Kutta, que são constituídos por uma família de métodos iterativos implícitos e explícitos. Os métodos de Runge-Kutta incluem o método de Euler, dada na equação (24), nesse caso temos o método de Runge-Kutta de primeira ordem e o método de Heun, equação (26), que é o método de Runge-Kutta de segunda ordem [²⁷]. Entre os métodos destacamos aqui o método de Runge-Kutta de quarta ordem (RK4), onde a vantagem desse método está em evitar os cálculos de derivadas de ordem superior, o qual exigiria um significativo esforço computacional além do que são problemas extremamentes complexos de se obter soluções analíticas.

No método RK4, a variável independente x da equação (22) é incrementada de δx (passo) em etapas e o novo valor da variável dependente y é calculado ao final de cada etapa de acordo com

(27)$$y_{n+1}=y_n+\frac{1}{6}\left(K_1+2K_2+2K_3+K_4\right)\delta x$$

e

(28)$$x_{n+1}=x_n+\delta x,$$

onde

(29)$$\begin{array}{lll}{K}_1\hfill & =\hfill & f\left(x_n,y_n\right),\hfill \\ K_2\hfill & =\hfill & f\left(x_n+\frac{1}{2}\delta x,y_n+\frac{1}{2}{K}_1\right),\hfill \\ K_3\hfill & =\hfill & f\left(x_n+\frac{1}{2}\delta x,y_n+\frac{1}{2}{K}_2\right)\text{e}\hfill \\ K_4\hfill & =\hfill & f\left(x_n+\delta x,y_n+K_3\right).\hfill \end{array}$$

Note que no conjunto das funções K em (29), K₂ depende de K₁ para determinar o valor de y no ponto x_n + δx/2 através do método de Euler, K₃ é dependente de K₂ e K₄ depende de K₃. Além disso, o conjunto de equações RK:Eq03 envolvem aproximações das derivadas no ponto médio do intervalo entre x_n e x_n+1.

A implementação computacional do método RK4 é apresentada em um fluxograma na Fig. 6 e é feita da seguinte forma: (i) entrar com as condições iniciais x = x₀ e y(x₀) = y₀, o passo δx e a função K₁ = f(x_n, y_n); (ii) é computado as funções K₂, K₃ e K₄ usando o método de Euler; (iii) em seguida é feita a estimativa de $\overline{f}(x_n,y_n)$; (iv) é calculado y_n+1; (v) enquanto x ˂ x_f as equações (27) e (28) são utilizadas para armazenar os dados; (vi) caso contrário, com os dados obtidos são plotados os gráficos ou guardados em uma lista (tabela de dados). O erro de truncamento local, erros causados pelo tipo de técnica empregada para a atualização do valor de y no determinado passo, do método RK4 é O(δx⁵), enquanto que para o método de Euler é O(δx²) e o método de Heun é de O(δx³).

Figura 6
Fluxograma método RK4.

4. Resultados e Discussão

Nessa seção apresentamos e discutimos, em duas partes, os resultados obtidos para os três modelos dos osciladores, Modelos 2.1, 2.2 e 2.3. A primeira parte, seção 4.1, trata-se do estudo numérico, onde é apresentado o comparativo entre os resultados provenientes de soluções exatas e soluções aproximadas obtidas através do métodos numéricos apresentados na seção 3. A segunda parte, seção 4.2, traz a discussão sobre a construção das animação para os simuladores visuais dos modelos.

4.1. Estudo numérico dos osciladores

Em todo estudo numérico dos osciladores, Modelos 2.1, 2.2 e 2.3, utilizamos três métodos numéricos implementado na linguagem de programação Python. O estudo foi desenvolvido na plataforma Anaconda Navigator e o Jupiter Lab versão 3.6.3, o qual utilizamos os pacotes NumPy, Pandas e Matplotilib como recursos fundamentais nesse trabalaho. A biblioteca NumPy é um pacote fundamental para computação científica em Python, ela fornece uma variedade de rotinas para operações rápidas de cálculos numéricos [³⁸]. Enquanto a biblioteca Pandas fornece estrutura de dados de alto desempenho e ferramentas de análise de dados, o qual ajuda a trabalahar com dados armazenados em planilhas ou bancos de dados [³⁹]. A bilioteca Matplotlib é utilizada para criação e visualizações de gráficos de dados em geral, essa biblioteca faz uso intenso de NumPy e outros códigos de extensão para fornecer bom desempenho [⁴⁰].

4.1.1. Estudo numérico do Modelo 2.1

Iniciamos o estudo numérico apartir do oscilador harmônico simples, Modelo 2.1, onde escrevemos a equação Eq02 na forma de duas equações diferenciais de primeira ordem, uma para a aceleração a(t) = dv(t)/dt = f₁(x, v, t) e outra para a velocidade v(t) = dx(t)/dt = f₂(x, v, t). Dessa forma, fazendo a mudança y → v e x → t temos que para o método de Euler (E) a equação (24) assume a seguinte forma

(30)$$v_{n+1}^E=v_n+f_1(x_n,v_n,t_n)\delta t$$

e a outra equação é obtida por y → x e x → t, também na equação ME:Eq03, assim

(31)$$x_{n+1}^E=x_n+f_2(x_n,v_n,t_n)\delta t.$$

As equações (30) e (31) são, respectivamente, as equações numéricas para velocidade e posição baseadas no método de Euler usadas no estudo do Modelo 2.1. Para o método de Heun (H), equação (25), as equações numéricas para o estudo são dadas por

(32)$$v_{n+1}^H=v_n+\frac{1}{2}\left[f_1(x_n,v_n,t_n)+f_1(x_{n+1},v_{n+1},t_{n+1})\right]\delta t,$$

e

(33)$$x_{n+1}^H=x_n+\frac{1}{2}\left[f_2(x_n,v_n,t_n)+f_2(x_{n+1},v_{n+1},t_{n+1})\right]\delta t.$$

Fazendo as mesmas mudanças no método RK4, temos as seguintes equações para a análise numérica

(34)$$v_{n+1}^{RK4}=v_n+\frac{1}{6}\left(K_{1v}+2K_{2v}+2K_{3v}+K_{4v}\right)\delta t$$

e

(35)$$x_{n+1}^{RK4}=x_n+\frac{1}{6}\left(K_{1x}+2K_{2x}+2K_{3x}+K_{4x}\right)\delta t.$$

Onde as funções K_iv, i = 1, 2, 3, e 4, são dadas por

$$\begin{array}{c}{K}_{1v}=f_1\left(x_n,v_n,t_n\right),\\ \begin{array}{ccc}{K}_{2v}& =& f_1\left(x_n+0.5K_{1x},v_n+0.5K_{1v},t_n+0.5\delta t\right),\\ K_{3v}& =& f_1\left(x_n+0.5K_{2x},v_n+0.5K_{2v},t_n+0.5\delta t\right),\end{array}\\ K_{4v}=f_1\left(x_n+K_{3x},v_n+K_{3v},t_n+\delta t\right).\end{array}$$

E para K_ix, i = 1, 2, 3 e 4, temos que

$$\begin{array}{c}{K}_{1x}=f_2\left(x_n,v_n,t_n\right),\\ \begin{array}{ccc}{K}_{2x}& =& f_2\left(x_n+0.5K_{1x},v_n+0.5K_{1v},t_n+0.5\delta t\right),\\ K_{3x}& =& f_2\left(x_n+0.5K_{2x},v_n+0.5K_{2v},t_n+0.5\delta t\right),\end{array}\\ K_{4x}=f_2\left(x_n+K_{3x},v_n+K_{3v},t_n+\delta t\right).\end{array}$$

Uma vez obtida as equações de posição e velocidade dos três métodos numéricos, equações (30) e (31) para o método de Euler, equações (32) e (33) para o método de Heun, e equações (34) e (35) para o método de RK4, usamos os parâmetros de massa m = 0.8, constante elástica da mola k = 4.0 e condições iniciais t = 0, x(0) = 2.0 e v(0) = 0.0 para obtemos os dados numéricos do modelo e compará-lo com o resultado analítico. Os passos de tempo δt = t_n+1 − t_n usados foram de 0.2 e 0.05 sendo respectivamente 50 e 200 interações para um parâmento de tempo máximo t = t_max = 10. O primeiro resultado que obtivemos é apresentado na Fig. 7, que mostra o plot das oscilações da posição x(t), Fig. 7(a), e da velocidade v(t), Fig. 7(b) para o passo δt = 0.2. Em ambas figuras, respectivamente posição x(t) e velocidade v(t), temos a solução analítica, equações (3) e (5), em preto (linha contínua), método de Euler, equações (32) e (33), em azul escuro (linha tracejada), método de Heun, equações Res:Eq04 e Res:Eq05, em vermelho (linha traço-ponto) e método de RK4, equações (34) e (35), em azul claro (linha pontilhada), para o passo δt = 0.2. Note que na Fig. 7 a solução numérica do método de RK4 praticamente coincide com o resultado analítico e uma melhor forma de observamos esses resultados é olhar a Tab. 2.

Tabela 2
Comparação para as posições x(t) do Modelo 2.1 obtidas analiticamente, equação (3), e pelos métodos de Euler, equação (31), Heun, equação (33), e RK4, equação (35), nos instantes de tempos t = 0, 3, 6 e 9. Na tabela temos, de cima para baixo, o passo δt = 0.2 e δt = 0.05 respectivamente.

t	Analítico	Euler	Heun	RK4	Erro % relativo Euler	Erro % relativo Heun	Erro % relativo RK4
δ⁢t = 0.2
0.00	2.000000	2.000000	2.000000	2.000000	0.000000	0.000000	0.000000
3.00	−0.649237	−0.900167	−0.540657	−0.649683	38.649985	16.724247	0.068696
6.00	−1.578491	−1.442437	−1.763005	−1.577655	8.619245	11.689265	0.052962
9.00	1.674053	1.793612	1.508791	1.674577	7.141889	9.871969	0.031301
δ⁢t = 0.05
0.00	2.000000	2.000000	2.000000	2.000000	0.000000	0.000000	0.000000
3.00	−0.649237	−0.716538	−0.642013	−0.649239	10.366168	1.112691	0.000308
6.00	−1.578491	−1.536118	−1.588673	−1.578488	2.684399	0.645046	0.000190
9.00	1.674053	1.711011	1.662236	1.674056	2.207696	0.705892	0.000179

Figura 7
Soluções analíticas e numéricas para o Modelo 2.1. Na Fig. (a) temos as soluções para a posição x(t)e em (b) as soluções para a velocidade v(t). Em ambas as figuras temos a solução analítica, equações (3) e (5), em preto (linha contínua), método de Euler, equações (30) e (31), em azul escuro (linha tracejada), método de Heun, equações (32) e (33), em vermelho (linha traço-ponto) e método de RK4, equações (34) e (35), em azul claro (linha pontilhada), para o passo δt = 0.2.

Embora utilizamos dois passos de tempo δt = 0.2 e δt = 0.05, no estudo numérico optamos por fazer oplot apenas para o δt = 0.2 pois os gráficos para os demais passos são semelhantes ao da Fig. 7 e do ponto de vista visual não é percebido a diferença dos resultados numéricos. Diante disso, construímos uma tabela de dados correspondentes aos valores da posição x(t), considerando seis casas decimais, obtidas pelos três métodos numéricos e pela solução analítica, como mostrado na Tab. 2. Perceba então que ao compararmos os dados obtemos que o método RK4 tem uma melhor aproximação numérica do resultado analítico.

Podemos ver isso na Tab. 2 escolhendo t = 6.00 para o passo δt = 0.2, por exemplo, onde o resultado analítico para a posição vale −1.578491,, enquanto que para o método de Euler temos −1.442437, para o método de Heun é 1.763005 e para RK4 temos 1.577655; e para o o passo δt = 0.05 temos a posição para o método de Euler dada por −1.536118,, para o método de Heun é −1.588673 e para RK4 temos −1.578488. A medida que diminuimos o passo a precisão dos resultados aumentam, mas isso requer mais cálculos computacionais.

Uma outra forma de verificarmos a precisão dos resultados além da comparação direta dos dados é analizar o erro percentual relativo, que é obtido calculando o módulo da diferença do valor real e o aproximado dividido pelo valor real multiplicado por 100, assim

(36)$$\text{Erro}\%\text{relativo =}\left|\frac{\text{valor real ࢤ aproximado}}{\text{valor real}}\right|\times \text{100},$$

o termo |valor real − aproximado| é chamado de valor absoluto. Temos então que o erro percentual relativo, equação (36), compara um valor real com o valor aproximado, expressando a resposta como o valor absoluto de uma porcentagem. Portanto, com a equação (36) foi possível calcular os erros relativos para os três métodos numéricos os quais obtemos que, na análise anterior em t = 6.00 e δt = 0.2, o Erro% relativo para o método de Euler é de 8.619245%, o Erro% relativo para o método de Heun é 11.689265% e o Erro% relativo para o método RK4 é de 0.052962% como mostrado na Tab. 2, respectivamente na sexta, sétima e oitava colunas. Isso mostra de fato que o método RK4 nos dá um resultado mais aproximado do resultado analítico quando comparado com os métodos de Euler e o de Heun. Na tentativa de diminuir os erros, diminuímos o passo para δt = 0.05, mesmo que isso quadruplicasse o número de cálculos necessários para garantirmos erros menores. Sendo assim, obtivemos para δt = 0.05 o Erro% relativo para o método de Euler foi de 2.684399%,, para o método de Heun foi de 0.645046% e para o método RK4 um erro de 0.000190%, esses resultados, para t = 6, discutidos aqui e outros, para t = 0,3 e t = 9, são apresentados na Tab. 2. Note então que, embora tenhamos um tempo de cálculo computacional maior, quando diminuímos os valores dos passos obtemos resultados ainda mais precisos nos três métodos utilizados.

4.1.2. Estudo numérico do Modelo 2.2

Prosseguindo com o nosso estudo numérico temos agora o modelo do pêndulo simples, Modelo 2.2, o qual escrevemos a equação mod2:Eq01 na forma de duas equações diferenciais de primeira ordem. Para a aceleração angular temos que α(t) = dω(t)/dt = f₁(θ, ω, t) e para a velocidade angular ω(t) = dθ(t)/dt = f₂(θ, ω, t) Dessa forma, fizemos a mudança y → ω e x → t tal que as equações numéricas baseadas no método de Euler usadas aqui são dadas por

(37)$${\omega }_{n+1}^E={\omega }_n+f_1({\theta }_n,{\omega }_n,t_n)\delta t,$$

e a outra equação é obtida por y → θ e x → t, também na equação (24), assim

(38)$${\theta }_{n+1}^E={\theta }_n+f_2({\theta }_n,{\omega }_n,t_n)\delta t.$$

Para o método de Heun, equação (25), as equações para o estudo são dadas por

(39)$${\omega }_{n+1}^H={\omega }_n+\frac{1}{2}\left[f_1({\theta }_n,{\omega }_n,t_n)+f_1({\theta }_{n+1},{\omega }_{n+1},t_{n+1})\right]\delta t,$$

e

(40)$${\theta }_{n+1}^H={\theta }_n+\frac{1}{2}\left[f_1({\theta }_n,{\omega }_n,t_n)+f_1({\theta }_{n+1},{\omega }_{n+1},t_{n+1})\right]\delta t.$$

Fazendo as mesmas mudança no método RK4, temos as seguintes equações para a análise numérica

(41)$${\omega }_{n+1}^{RK4}=v_n+\frac{1}{6}\left(K_{1\omega }+2K_{2\omega }+2K_{3\omega }+K_{4\omega }\right)\delta t$$

e

(42)$${\theta }_{n+1}^{RK4}={\theta }_n+\frac{1}{6}\left(K_{1\theta }+2K_{2\theta }+2K_{3\theta }+K_{4\theta }\right)\delta t.$$

Com as funções K_iω, i = 1, 2, 3 e 4, dadas por

(57)$$\begin{array}{c}{K}_{1\omega }=f_1\left({\theta }_n,{\omega }_n,t_n\right),\\ \begin{array}{ccc}{K}_{2\omega }& =& f_1\left({\theta }_n+\frac{1}{2}{K}_{1\theta },{\omega }_n+\frac{1}{2}{K}_{1\omega },t_n+\frac{1}{2}\delta t\right),\\ K_{3\omega }& =& f_1\left({\theta }_n+\frac{1}{2}{K}_{2\theta },{\omega }_n+\frac{1}{2}{K}_{2\omega },t_n+\frac{1}{2}\delta t\right),\end{array}\\ K_{4\omega }=f_1\left({\theta }_n+K_{3\theta },{\omega }_n+K_{3\omega },t_n+\delta t\right).\end{array}$$

E para K_iθ, i = 1, 2, 3 e 4, temos que

(58)$$\begin{array}{c}{K}_{1\theta }=f_2\left({\theta }_n,{\omega }_n,t_n\right),\\ \begin{array}{ccc}{K}_{2\theta }& =& f_2\left({\theta }_n+\frac{1}{2}{K}_{1\theta },{\omega }_n+\frac{1}{2}{K}_{1\omega },t_n+\frac{1}{2}\delta t\right),\\ K_{3\theta }& =& f_2\left({\theta }_n+\frac{1}{2}{K}_{2\theta },{\omega }_n+\frac{1}{2}{K}_{2\omega },t_n+\frac{1}{2}\delta t\right),\end{array}\\ K_{4\theta }=f_2\left({\theta }_n+K_{3\theta },{\omega }_n+K_{3\omega },t_n+\delta t\right).\end{array}$$

Com as equações de posição e velocidade dos três métodos numéricos, equações Res:(37) e (38) para o método de Euler, equações (39) e (40) para o método de Heun, e equações (41) e (42) para o método de RK4, usamos, em unidades de medidas abritárias, os parâmetros de massa da barra m_b = 1.0, massa do disco m_d = m_b, raio do disco R = 6, comprimento da barra L = 4.0R, a dis tância da mola ao pivô d = L/3, a contantes elástica da mola k = 6.0, a aceleração da gravidade g = 9.8 e as condições iniciais t = 0, θ(0) = π/18 e ω(0) = 0.0 para obtemos os dados numéricos do modelo e compará-lo com o resultado analítico. Assim como feito no Modelo 2.1, o passo de tempo δt = t_n+1 − t_n usado foi de 0.2 e 0.05 sendo respectivamente 50 e 200 interações para um parâmento de tempo máximo t = t_max = 10. Com esses parâmetros plotamos a Fig. 8 que mostra as oscilações da posição angular θ(t), Fig. 8(a), e da velocidade angular ω(t), Fig. 8(b) para o passo δt = 0.2. Em ambas as figuras temos a solução analítica, equações (11) e (13), em preto (linha contínua), método de Euler, equações (37) e (38), em azul escuro (linha tracejada), método de Heun, equações (39) e (40), em vermelho (linha traço-ponto) e método RK4, equações (41) e (42), em azul claro (linha pontilhada), para o passo δt = 0.2. Assim como o Modelo 2.1, estudado anteriomente, tivemos aqui que a solução numérica para o método RK4 apresentou resultados bem mais próximos do resultado analítico quando comparado com os métodos de Euler e o de Heun, ver Fig. 8.

Figura 8
Soluções analíticas e numéricas para o Modelo 2.2. Na Fig. (a) temos as soluções para a posição x(t) e em (b) as soluções para a velocidade v(t). Em ambas as figuras temos a solução analítica, equações (11) e (13), em preto (linha contínua), método de Euler, equações (37) e (38), em azul escuro (linha tracejada), método de Heun, equações (39) e (40), em vermelho (linha traço-ponto) e método de RK4, equações (41) e (42), em azul claro (linha pontilhada), para o passo δt = 0.2.

Para verificarmos detalhadamente os resultados numéricos construímos também, assim como no Modelo 2.1, uma tabela de dados correspondentes aos valores da posição angular θ(t) para três passos de tempo δt = 0.2 e δt = 0.05 como mostra a Tab. 3. Como a quantidade de dados é grande, escolhemos o instante de tempo t = 6.00 para discutimos aqui a comparação dos três métodos numéricos nos três passos de tempo, para os outros valores a discussão é similar. Perceba então que para o passo δt = 0.2 onde o resultado analítico para a posição angular vale 0.05114, enquanto que para o método de Euler temos 0.026461, para o método de Heun é 0.038498 e para RK4 temos 0.051179; e para o o passo δt = 0.05 temos a posição para o método de Euler dada por 0.045673, para o método de Heun é 0.050319 e para RK4 temos 0.05114. Note que já no passo δt = 0.2 o método RK4 já apresenta um resultado bem próximo a solução analítica e isso fica ainda mais evidente quando analisamos o erro percentual relativo pela equação (36). Sendo assim, obtemos que em t = 6.00 e δt = 0.2, o Erro% relativo para o método de Euler é de 48.257724%, o Erro% relativo para o método de Heun é 24.720375% e o Erro% relativo para o método RK4 é de 0.076261%; E por fim para δt = 0.05 o Erro% relativo para o método de Euler foi de 10.690262%, para o método de Heun foi de 1.605397% e para o método RK4 um erro de 0.000000%, como mostrado na Tab. 3, respectivamente na sexta, sétima e oitava colunas. Todos os resultados, para t = 4, discutidos aqui e outros, para t = 0, 3, 6 e 9, são apresentados na Tab. 3 o qual mostra que o método RK4 fornece uma melhor aproximação do resultado analítico em todos os passos quando comparado com os métodos de Euler e Heun.

Tabela 3
Comparação para as posições θ do Modelo 2.2 obtidas analiticamente, equação (11), e pelos métodos de Euler, equação (37), Heun, equação (39), e RK4, equação (41), nos instantes de tempos t = 0, 3, 6 e 9. Na tabela temos, de cima para baixo, o passo δt = 0.2 e δt = 0.05 respectivamente.

t	Analítico	Euler	Heun	RK4	Erro % relativo Euler	Erro % relativo Heun	Erro % relativo RK4
δ⁢t = 0.2
0.00	0.174533	0.174533	0.174533	0.174533	0.000000	0.000000	0.000000
3.00	−0.140334	−0.125931	−0.137193	−0.140343	10.263372	2.238232	0.006413
6.00	0.05114	0.026461	0.038498	0.051179	48.257724	24.720375	0.076261
9.00	0.058095	0.083697	0.078755	0.058029	44.069197	35.562441	0.113607
δ⁢t = 0.05
0.00	0.174533	0.174533	0.174533	0.174533	0.000000	0.000000	0.000000
3.00	−0.140334	−0.136998	−0.140091	−0.140334	2.377186	0.173158	0.000000
6.00	0.05114	0.045673	0.050319	0.05114	10.690262	1.605397	0.000000
9.00	0.058095	0.063585	0.059346	0.058095	9.450039	2.153369	0.000000

4.1.3. Estudo numérico do Modelo 2.3

A construção da análise numérica do oscilador acoplado, Modelo 2.3, é mais sofistificada que os Modelos 2.1 e 2.2 pois aqui temos um par de equações diferencais de segunda ordem acopladas, equações mod3:Eq01 e mod3:Eq02. Sendo assim, escrevemos essas equações na forma de quatro equações diferenciais de primeira ordem dadas por: aceleração e velocidade do corpo 1, respectivamente, a₁(t) = dv₁(t)/dt = f₁₁(x₁, x₂, v₁, v₂, t), v₁(t) = dx₁(t)/dt = f₂₁(x₁, x₂, v₁, v₂, t) e aceleração e velocidade para o corpo 2, temos respectivamente a₂(t) = dv₂(t)/dt = f₁₂(x₁, x₂, v₁, v₂, t), v₂(t) = dx₂(t)/dt = f₂₂(x₁, x₂, v₁, v₂, t). Então, para o estudo desse modelo escrevemos as equações numéricas baseadas no método de Euler, equação (24), nas formas

(43)$$v_{1_{n+1}}^E=v_{1_n}+f_{11}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\delta t,$$

(44)$$x_{1_{n+1}}^E=x_{1_n}+f_{21}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\delta t,$$

(45)$$v_{2_{n+1}}^E=v_{2_n}+f_{12}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\delta t,$$

(46) $$x_{2_{n+1}}^E=x_{2_n}+f_{22}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\delta t,$$

essas equações foram obtidas fazendo as seguintes alterações no método de Euler: para a velocidade y → v₁ e y → v₂; para a posição y → x₁ e y → x₂; e para o tempo x → t. Para o método de Heun, equação (25), as equações numéricas para o estudo são dadas por

(47)$$v_{1_{n+1}}^H=v_{1_n}+\frac{1}{2}\left[f_{11}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\right.+\left.f_{11}(x_{1_{n+1}},x_{2_{n+1}},v_{1_{n+1}},v_{2_{n+1}},t)\right]\delta t,$$

(48)$$x_{1_{n+1}}^H=x_{1_n}+\frac{1}{2}\left[f_{21}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\right.+\left.f_{11}(x_{1_{n+1}},x_{2_{n+1}},v_{1_{n+1}},v_{2_{n+1}},t)\right]\delta t,$$

(49)$$v_{2_{n+1}}^H=v_{2_n}+\frac{1}{2}\left[f_{12}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\right.+\left.f_{11}(x_{1_{n+1}},x_{2_{n+1}},v_{1_{n+1}},v_{2_{n+1}},t)\right]\delta t,$$

(50)$$x_{2_{n+1}}^H=x_{2_n}+\frac{1}{2}\left[f_{22}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t)\right.+\left.f_{11}(x_{1_{n+1}},x_{2_{n+1}},v_{1_{n+1}},v_{2_{n+1}},t)\right]\delta t,$$

Fazendo as mesmas mudança no método RK4 dado pelas equações (27) e (29), temos as seguintes equações para a análise numérica

(51)$$v_{1_{n+1}}^{RK4}=v_{1_n}+\frac{1}{6}\left(K_{11v}+2K_{21v}+2K_{31v}+K_{41v}\right)\delta t,$$

(52)$$x_{1_{n+1}}^{RK4}=x_{1_n}+\frac{1}{6}\left(K_{11x}+2K_{21x}+2K_{31x}+K_{41x}\right)\delta t,$$

(53)$$v_{2_{n+1}}^{RK4}=v_{2_n}+\frac{1}{6}\left(K_{12v}+2K_{22v}+2K_{32v}+K_{42v}\right)\delta t,$$

(54)$$x_{2_{n+1}}^{RK4}=x_{2_n}+\frac{1}{6}\left(K_{12x}+2K_{22x}+2K_{32x}+K_{42x}\right)\delta t.$$

As funções K_ijv, com i = 1, 2, 3 e 4; e j = 1 e 2, são dadas por:

(59)$$\begin{array}{lll}{K}_{11v}\hfill & =\hfill & f_{11}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t),\hfill \\ K_{21v}\hfill & =\hfill & f_{11}\left(x_{1_n}+\frac{1}{2}{K}_{11x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}+\left.\frac{1}{2}{K}_{11v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\right.\hfill \\ K_{31v}\hfill & =\hfill & f_{11}\left(x_{1_n}+\frac{1}{2}{K}_{21x},x_{2_n}+\frac{1}{2}{K}_{31x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{21v},v_{2_n}+\frac{1}{2}{K}_{31v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{41v}\hfill & =\hfill & f_{11}\left(x_{1_n}+K_{31x},x_{2_n}+K_{32x},v_{1_n}\right.+\left.K_{31v},v_{2_n}+K_{32v},t_n+\delta t\right);\hfill \\ K_{12v}\hfill & =\hfill & f_{12}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t),\hfill \\ K_{22v}\hfill & =\hfill & f_{12}\left(x_{1_n}+\frac{1}{2}{K}_{11x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{11v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{32v}\hfill & =\hfill & f_{12}\left(x_{1_n}+\frac{1}{2}{K}_{21x},x_{2_n}+\frac{1}{2}{K}_{31x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{21v},v_{2_n}+\frac{1}{2}{K}_{31v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{42v}\hfill & =\hfill & f_{12}\left(x_{1_n}+K_{31x},x_{2_n}+K_{32x},v_{1_n}\right.+\left.K_{31v},v_{2_n}+K_{32v},t_n+\delta t\right);\hfill \end{array}$$

E para as funções K_ijx, com i = 1, 2, 3 e 4; e j = 1 e 2, temos: e 2, temos:

(60)$$\begin{array}{lll}{K}_{11x}\hfill & =\hfill & f_{21}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t),\hfill \\ K_{21x}\hfill & =\hfill & f_{21}\left(x_{1_n}+\frac{1}{2}{K}_{11x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{11v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{31x}\hfill & =\hfill & f_{21}\left(x_{1_n}+\frac{1}{2}{K}_{21x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{21v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{41x}\hfill & =\hfill & f_{22}\left(x_{1_n}+K_{31x},x_{2_n}+K_{32x},v_{1_n}\right.+\left.K_{31v},v_{2_n}+K_{32v},t_n+\delta t\right)\hfill \\ K_{12x}\hfill & =\hfill & f_{22}(x_{1_n},x_{2_n},v_{1_n},v_{2_n},t),\hfill \\ K_{22x}\hfill & =\hfill & f_{22}\left(x_{1_n}+\frac{1}{2}{K}_{11x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{11v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{32x}\hfill & =\hfill & f_{22}\left(x_{1_n}+\frac{1}{2}{K}_{21x},x_{2_n}+\frac{1}{2}{K}_{21x},v_{1_n}\right.+\left.\frac{1}{2}{K}_{21v},v_{2_n}+\frac{1}{2}{K}_{21v},t_n+\frac{1}{2}\delta t\right),\hfill \\ K_{42x}\hfill & =\hfill & f_{22}\left(x_{1_n}+K_{31x},x_{2_n}+K_{32x},v_{1_n}\right.+\left.K_{31v},v_{2_n}+K_{32v},t_n+\delta t\right)\hfill \end{array}$$

No estudo, consideramos por simplicidade o caso de molas iguais, k₁ = k₂ = k, então as frequêcias ω₁ e ω₂, dadas respectivamente pelas equações 1mod3:Eq05 e 2mod3:Eq05, assumiram a forma ${\omega }_1={\left(3-\sqrt{5}\right)}^{1/2}\sqrt{k/2m}$ e ${\omega }_2={\left(3+\sqrt{5}\right)}^{1/2}\sqrt{k/2m}$. Consideramos também que os coeficientes da amplitude de movimento A₁ = A₂ = 1, sendo assim os coeficientes B₁ e B₂ são dados por $B_1=\left(1+\sqrt{5}\right)/2$ e $B_2=\left(1-\sqrt{5}\right)/2$. Os parâmetros numéricos usados foram: massa m = 0.8, constante elástica da mola k = 1.5 e condições iniciais para o corpo 1 x₁(0) = A₁ + A₂ = 2.0 e v₁(0) = 0.0 e para o corpo 2 x₂(0) = B₁ + B₂ = 1.0 e v₂(0) = 0.0. O passo de tempo δt = t_n+1 − t_n usado foi de 0.2 e 0.05 sendo respectivamente 75 e 300 interações para um parâmento de tempo máximo t = t_max = 15.

Com os parâmetros e condições inicias dadas, plotamos na Fig. 9 as soluções analíticas e numéricas que mostram as oscilações para as posições e velocidades dos corpos 1 e 2 para o passo δt = 0.05. Na Fig. 9(a) temos as soluções para a posição x₁(t), figura de cima, e velocidade v₁(t), figura de baixo, para o corpo 1 apresentadas na forma: solução analítica, equações (16) e (20), em preto (linha contínua); método de Euler, equações (44) e (43), em azul escuro (linha tracejada); método de Heun, equações (48) e (47), em vermelho (linha traço-ponto); e método de RK4, equações (52) e (51), em azul claro (linha pontilhada). A ordem das equações são respectivamente dadas para a posição x₁ e a velocidade v₁. A Fig. 9(b) apresenta as soluções para a posição x₂(t), figura de cima, e velocidade v₂(t), figura de baixo, para o corpo 2 na seguinte forma: solução analítica, equações (16) e (21), em preto (linha contínua); método de Euler, equações. (46) e (45), em azul escuro (linha tracejada); método de Heun, equações (50) e (49), em vermelho (linha traço-ponto); e método de RK4, equações (54) e (53), em azul claro (linha pontilhada). A ordem das equações são respectivamente dadas para a posição x₁ e a velocidade v₁.

Figura 9
Soluções analíticas e numéricas para o Modelo 2.3. Na Fig. (a) temos as soluções para a posição x₁(t) e velocidade v₁(t) para o corpo 1, sendo a solução analítica, equações (16) e (20), em preto (linha contínua), método de Euler, equações (44) e (43) respectivamente, em azul escuro (linha tracejada), método de Heun, na ordem equações (48) e (47), em vermelho (linha traço-ponto) e método de RK4, respectivamente dadas por equações (52) e (51), em azul claro (linha pontilhada). A Fig (b) apresenta as soluções para a posição x₂(t)e velocidade v₂(t) para o corpo 2, sendo a solução analítica, equações (17) e (21) respectivamente, em preto (linha contínua), método de Euler, na ordem equações (46) e (45), em azul escuro (linha tracejada), método de Heun, equações (50) e (49), em vermelho (linha traço-ponto) e método de RK4, respectivamente dadas pelas equações (54) e (53), em azul claro (linha pontilhada). Em ambas figuras temos o passo δt = 0.05.

Perceba na Fig. 9 que visualmente os três métodos apresentam boas aproximações até t = 5, acima disso os resultados numéricos começam a divergir significativamente, isso é bem notável para passos de tempo maiores que δt = 0.05. Sendo assim, construímos uma tabela de dados correspondentes aos valores das posições x₁ e x₂ para analisarmos os comportamentos dos métodos numéricos para os passos δt = 0.2 e δt = 0.05, como mostrado na Tab. 4. Na discussão a seguir, escolhemos o instante de tempo t = 12 para a comparação dos três métodos numéricos nos três passos de tempo. Para os demais valores a discussão é similar e está também apresentada na Tab. 4.

Tabela 4
Comparação para as posições x₁ e x₂ do Modelo 2.3 obtidas analiticamente, equações (16) e (17), e pelos métodos de Euler, equação (37), Heun, equação (39), e RK4, equação (41), nos instantes de tempos t = 0, 2, 4, 6, 8 e 10. Na tabela temos, de cima para baixo, os passos de tempo δt = 0.2 e δt = 0.05 para as posições x₁ e x₂ respectivamente.

t	Analítico	Euler	Heun	RK4	Erro % relativo Euler	Erro % relativo Heun	Erro % relativo RK4
x1 para δ⁢t = 0.2
0.00	2.000000	2.000000	2.000000	2.000000	0.000000	0.000000	0.000000
3.00	0.110859	−0.636435	0.140181	0.127227	674.094120	26.449815	14.764701
5.00	−0.380411	−0.018542	0.253902	−0.123560	95.125798	166.744127	67.519341
7.00	−0.044126	1.360077	−0.624686	−0.307866	3182.257626	1315.686897	597.697503
10.00	−1.558350	−2.491956	−2.707683	−2.915612	59.909905	73.753201	87.096095
12.00	−0.628540	−1.307307	−2.829227	−1.169905	107.991059	350.126802	86.130557
15.00	0.747309	3.174851	−2.615211	−0.300278	324.837785	449.950422	140.181237
x1 para δ⁢t = 0.05
0.00	2.000000	2.000000	2.000000	2.000000	0.000000	0.000000	0.000000
3.00	0.110859	−0.065405	0.107091	0.107241	158.998367	3.398912	3.263605
5.00	−0.380411	−0.319038	−0.301570	−0.318872	16.133340	20.725216	16.176977
7.00	−0.044126	0.270351	−0.114825	−0.105557	712.679599	160.220732	139.217242
10.00	−1.558350	−1.598460	−1.808284	−1.811058	2.573876	16.038374	16.216383
12.00	−0.628540	−0.807988	−0.786556	−0.739862	28.549973	25.140166	17.711204
15.00	0.747309	1.139183	0.563394	0.623195	52.438014	24.610302	16.608123
x2 para δ⁢t = 0.2
0.00	1.000000	1.000000	1.000000	1.000000	0.000000	0.000000	0.000000
3.00	−1.910531	−2.521425	−2.327695	−2.341878	31.975090	21.834977	22.577336
5.00	−0.799444	−0.352885	−0.936559	−0.698019	55.858697	17.151295	12.686942
7.00	2.120598	3.446956	3.262449	3.076594	62.546414	53.845708	45.081435
10.00	−0.315652	−3.262170	−0.254298	−0.226934	933.470404	19.437228	28.106269
12.00	−1.276971	−1.799743	−0.302942	−1.340084	40.938439	76.276517	4.942399
15.00	1.755971	5.321681	4.581296	2.839654	203.062010	160.898158	61.714174
x2 para δ⁢t = 0.05
0.00	1.000000	1.000000	1.000000	1.000000	0.000000	0.000000	0.000000
3.00	−1.910531	−2.039323	−2.001033	−2.003910	6.741163	4.737008	4.887594
5.00	−0.799444	−0.764301	−0.800616	−0.789802	4.395930	0.146602	1.206088
7.00	2.120598	2.354112	2.322067	2.317896	11.011705	9.500575	9.303885
10.00	−0.315652	−0.833255	−0.346837	−0.344717	163.979002	9.879551	9.207925
12.00	−1.276971	−1.399589	−1.252715	−1.280279	9.602254	1.899495	0.259051
15.00	1.755971	2.322377	2.036675	1.996691	32.256000	15.985685	13.708655

Em relação a posição x₁ em t = 12, para o passo δt = 0.2, temos o resultado analítico dado por −0.628540, enquanto que para o método de Euler −1.307307, para o método de Heun é −2.829227 e para RK4 temos −1.169905. Note então que esses resultados não são boas aproximações e isso é comprovado pelo o erro percentual relativo calculado na equação erro1 o qual mostra os erros% para os métodos de Euler, Heun e RK4, respectivamente, dados por 107.991059%, 350.126802% e 86.130557%, como apresentado na sexta, sétima e oitava coluna da Tab. 4. Para um passo de tempo ainda menor, δt = 0.05, obtemos resultados mais próximos do analítico para os três métodos, são eles: −0.807988, −0.786556 e −0.739862, nessa ordem, para os métodos de Euler, Heun e RK4, enquando que os erros% relativo para os métodos de Euler, Heun e RK4 são dados respectivamente por 28.549973%, 25.140166% e 17.711204%.

Analisando agora a posição x₂ em t = 12, para o passo δt = 0.2, obtemos o resultado analítico dado por −1.276971, para o método de Euler −1.799743 o que nos dá 40.938439% para o erro% relativo, no método de Heun temos que a posição é −0.302942 e para RK4 temos −1.340084 para posição e erro% relativo de 4.942399%. Por último temos que para δt = 0.05 obtemos os seguintes resultados numéricos: −1.399589, −1.252715 e −1.280279, respectivamente, para os métodos de Euler, Heun e RK4, enquando que os erros% relativo para o método de Euler foi de 9.602254%, método de Heun 1.899495%> e método RK4 0.259051%. Perceba que, nos resultados obtidos, temos o método RK4 como sendo o método numérico que mais nos fornece uma boa aproximação do resultado analítico como mostra a Tab. 4.

Vale ressaltar também, que como estamos diante de um sistema acoplado é natural obtermos resultados numéricos que divergem significativamente dos resultados analíticos, como ocorreu para x₁, ver Tab. 4, no passo δt = 0.2, pois temos acúmulos de erros sucessivos na aplicação dos métodos numéricos de uma equação diferencial para outra. Diante disso, é preciso diminuir cada vez mais o passo de tempo δt para obtermos resultados numéricos mais aproximados dos resultaos analíticos. Outro ponto importante de mencionar, é que a escolha de apresentar os resultados numéricos para a posição nas Tabelas 2, 3 e 4 nos três modelos discutidos acima se dá pelo fato de que o comportamento oscilatorio da velocidade é similar ao da posição a menos de uma constante.

4.2. Construção e visualização dos osciladores

Na construção e visualização da animação dos osciladores, Modelos 2.1, 2.2 e 2.3, utilizamos linguagem de programação Python 3.12.4 e a biblioteca VPython 3.2 − 5.74, ambos instaladores disponíveis respectivamente nas Refs. [⁴¹, ⁴²]. A biblioteca VPython é um ambiente poderoso e fácil de usar para criar animações 3D navegáveis em tempo real através de um sistema de coordenadas x, y, z. Um ponto interessante dos recursos que o VPython oferece é a utilização dos próprios objetos como variáveis, o qual nos proporciona a criar e utilizar uma infinidade de objetos e formas integrados, bem como recursos gráficos, adicionando propriedades físicas como tempo, posição e massa, o que torna o programa mais interativo [⁴³, ⁴⁴].

Em nosso estudo para os Modelos 2.1, 2.2 e 2.3, iniciamos primeiro o processo para tornar o programa interativo começando com a criação da cena, objetos e visualização dos gráficos. Para isso, é necessário importar o pacote visual que é responsável pelos recursos de programação para a visualização 3D, então inserimos o comando from visual import * (aparece com essas cores na tela do programa), onde, o símbolo * indica que deseja-se ter acesso a todos os recursos do pacote. Inserimos também os comandos from visual.graph import * e from math import * que são os pacotes de visualização dos gráficos e de acesso às funções matemática, respectivamente.

Para a criação da cena, dimensões da tela de visualização da animação, usamos as seguintes funções: scene.range (zoom em relação ao centro da tela), scene.center (localização do centro da tela), scene.background (cor de fundo da tela), scene.width (largura da cena) e scene.height (altura da cena). Em seguida, construímos os objetos para cada um dos modelos usando várias funções que dão as formas geométricas de cada objetos, essas funções serão apresentadas em cada um dos modelos nos parágrafos a seguir.

Em relação a visualização dos gráficos usamos a função gdots() o qual exibe na tela uma lista de pares (x, y) como pontos discretos plotando então a curva de interesse no estudo. Associamos uma cena com a função gdots() para apresentar os gráficos, para isso criamos uma cena separada com a função gdisplay() com os seguinte comandos: largura (width), altura (height), título na tela (title), nomes nas coordenadas x (xtitle) e y (ytitle) além de um itervalo de apresentação de dados para cada uma das coordenadas (xmax, xmin) e (ymax, ymin) e cor de fundo (background=color.white), nesse caso está uma cor branca. Por fim, combinamos as funções gdots() e gdisplay() na forma gdots(gdisplay=grafico). Por exemplo, podemos contruir uma cena da seguinte forma

Neste exemplo de cena, a janela do gráfico será localizado em (0,0), com um tamanho de 800 por 400 pixels, e na barra de título vai constar Posição vs. Tempo. O gráfico terá um título de t(s) no eixo horizontal e x(m) no eixo vertical, terá limites fixados de 0 até 10 no eixo horizontal se estenderá de −3 para 3 no eixo vertical. A cor de primeiro plano (branco por padrão) é preta, e a cor de fundo (preto por padrão) é branca. Na cena está presente um gráfico de cor azul, como mostrado na Fig. 10.

Figura 10
>Exemplo de cena construída utilizando algumas funções da biblioteca VPython.

A construção e tradução dos modelos físicos, Modelos 2.1, 2.2 e 2.3, para a linguagem computacional foi feita usando a segunda Lei de Newton $a=\sum F/m$ combinada com o método de Euler, equação (24), para as equações de velocidade e posição dadas respectivamente por v_n+1 = v_n + aδt e x_n+1 = x_n + v_nδt, onde δt é o passo de tempo e a força resultante ($\sum F$) que atua no sistema vai depender de cada um dos modelos estudados. Uma vantagem aqui na tradução do sistema físico para a linguagem computacional é o fato de que não precisamos nos preocupar com a notação de equação diferencial, isso facilita muito na escrita e no entendimneto dos programas para quem não é familiarizado com tal notação.

Para concluir a animação dos modelos, introduzimos a função while() que é uma instrução de fluxo de controle permitindo que o código seja executado repetidamente, dependendo se uma determinada condição é satisfeita ou não. Quando a condição se torna falsa, a linha imediatamente após o loop do programa é executada. Usamos também o comando rate() que é uma parte fundamental de qualquer loop de animação em um programa VPython, pois com ele controlamos a velocidade com que são feitas as iterações. O rate() permite que o VPython redesenhe a exibição, para que os objetos possam se mover, as animações possam ocorrer e os gráficos possam ser plotados. Sem esse comando rate() o programa seria executado mais rapidamente [⁴⁵]. Por fim, incluimos a função condicional if que é uma das instruções condicionais mais usadas em linguagens de programação. Essa condição decide se certas instruções precisam ser executadas ou não. Sendo assim, é verificado em uma determinada condição, se a condição for verdadeira, então o conjunto de códigos presente dentro do bloco if será executado, caso contrário não.

4.2.1. Construção e visualização do Modelo 2.1

Uma vez apresentada as funções e como foi feita a tradução dos modelos físicos para a linguagem computacional vamos agora mostrar a construção para cada modelo separado. No sistema massa-mola com rotação, Modelo 2.1, construímos duas cenas: cena 1 - visualização da dinâmica do oscilador; e cena 2 - apresentação dos gráficos da posição x(t), velocidade v(t) aceleração a(t). Para a contrução da cena 1 definimos os objetos e a atribuição de variáveis a eles. Sendo assim, usamos a função box() para desenharmos uma parede vertical esquerda, um piso e um suporte. Para a função box() usamos os comandos pos=(), size=() e color=color.– que são responsáveis por definirem a posição, tamanho e cor dos objetos. Utilizamos a função cylinder() para desenharmos um disco e um eixo que passa pelo centro do cilindro, para essa função usamos os comandos pos=(), radius=()axis=() e color=color.– que definem a posição, o raio, o eixo e cor dos objetos. Para completar os objetos desse modelo, desenhamos também uma mola inserindo a função helix() em conjunto com os comandos pos=vector(-7,0,0), axis=vector(-1, 4, 0), radius=0.5, color=color.orange, thickness=0.1, coils=10) e a função arrow=() para desenharmos uma seta que indica o sentido da rotação do cilindro na animação.

Em relação ao sistema físico-computacional consideramos os seguintes parâmetros em unidades arbitrárias: massa m = 0.4,, constante da mola k = 4.0, raio do cilindro R = 0.1; condições iniciais x(0) = 2.0 e v(0) = 0.0; e passo de tempo δt = 0.01. Com esses parâmetros geramos a animação do Modelo 2.1 usando a função while() e escrevendo as equações de movimento para a aceleração a = .(2k/3m) . disco.pos.x, como mostrado na equação Eq02, velocidade v + = a ∗ dt, posição x + = v∗dt e velocidade angular ω + = a∗dt/R, onde o termo disco.pos.x é a posição x do disco e o operador + = é adicionada para somar o próximo valor a um valor já existente, ou seja, é um acumulador de dados, então x + = n é uma abreviação de x = x + n

A animação gerada para a mola se dá através de mola.axis = suporte3.pos − paredeE.pos + (0, 1, 0), essa expressão descreve o movimento da mola que está entre as consecutivas posições do suporte preso ao disco e a posição da parede fixa. Para a rotação do disco é obtida usando a expressão disco.rotate(axis = (0, 0, 0.5), angle = −w ∗ dt), onde descreve o rolamento do disco na superfície plana.

Algumas posições específicas para animação (conjunto: disco, mola, parede vertical esquerda, suporte e piso) podem ser vistas em x = 2, Fig. 11(a), x = 0, Fig. 11(b), e x = −2, Fig. 11(c). Para finalizar a construção desse modelo usamos a função condicional if para limitarmos a visualização da animação dos gráficos da posição x(t), velocidade v(t) e aceleração a(t), respectivamente nas cores vermelho, preto e zul, como mostrado na cena 2, Fig. 11(d) para t ˂ 10, onde a evolução do parâmentro de tempo foi calculada por t + = dt. O código completo utilizado na construção do Modelo 2.1 pode ser visto no Apêndice [codigo1]B.

Figura 11
Visualização da contrução do Modelo 2.1 para as posições x = 2, Fig. (a), x = 0, Fig. (b), e x = −2, Fig. (c), onde todas essas posições estão presentes na cena 1 nomeada de Sistema Massa-mola com Rotação. Temos também os gráficos do movimento para a posição x(t), velocidade v(t) e aceleração a(t), respectivamente nas cores vermelho, preto e zul, apresentado na Fig. (d) para a cena 2 nomeada de Gráficos para o movimento.

4.2.2. Construção e visualização do Modelo 2.2

A construção da animação para o pêndulo físico acoplado, Modelo 2.2, foi feita através de duas cenas: cena 1 – visualização da dinâmica do modelo, nomeada de Pêndulo físico acoplado; e cena 2 – apresentação dos gráficos da posição angular θ(t), velocidade angular ω(t)e aceleração angular α(t), nomeada de Gráficos para o movimento.

Na a contrução da cena 1, definimos os objetos a partir de funções já apresentadas no Modelo 2.1. Sendo assim, usamos a função box() para desenharmos a parede esquerda, o teto e a barra vertical do pêndulo físico junto ao comando material=materials.wood que é responsável por simular uma textura no objeto de interesse, aqui usamos uma textura de madeira como mostrado na Fig. 12. Adicionamos a função cylinder() e os comandos pos=(3.5,4.4,−1.5), radius=0.1, axis=(0,0,-2) e color = color.red, que indicam respectivamente a posição, o raio, eixo e a cor, para criarmos um pino de fixação de referência para o ponto de rotação. Também usamos a mesma função cylinder() para criarmos um disco fixado a outra extremidade da barra. Por fim, desenhamos uma mola inserindo a função helix() com o comando material=materials.rough no programa.

A animação na cena 1 foi contruída através da tradução do sistema físico-computacional. Para isso, consideramos os seguintes parâmetros em unidades arbitrárias: massas m_d = m_b = 1.0,, raio do disco R = 6, comprimento da barra L = 4R, distância vertical da mola ao pivô d = L/3, constante elástica da mola k = 6 e aceleração da gravidade g = 9.8. Para as condições iniciais consideramos que a posição angular inicial theta(0) = π/18 e velocidade angular inicial omega(0) = 0.0, com o passo de tempo δt = 0.01. Usamos a função while() e escrevemos as equações de movimento, considerando pequenas oscilações, para a aceleração angular alpha = −(w₀ ∗ w₀) ∗ theta, onde w₀ = sqrt((m_b ∗ g ∗ L + 2 ∗ md ∗ g ∗ (L + R) + 2 ∗ k ∗ d ∗ d)/(2∗I)), velocidade angular omega + = alpha∗dt e posição angular theta + = omega ∗ dt.

O movimento do pêndulo é gerado usando que x = 4.5 ∗ theta, barra.pos.x = disco.pos.x, mola.axis = barra.pos−paredeE.pos+(−3.5,−1, 0), disco.pos.x = x e pendulo.rotate(axis = (0, 0, 1), angle = omega ∗ dt) que descrevem, respectivamente, a posição angular, o movimento da barra, o movimento da mola, a oscilação do disco e a rotação de todo o pêndulo em torno do ponto de pivô. A Fig. 12 apresenta algumas posições específicas para animação do Modelo 2.2, onde temos a condição inicial θ = π/18, Fig. 12(a), temos a posição θ = −π/18, Fig. 12(b), e por fim a posição mais baixa da trajetória θ = Fig. 12(c).

Figura 12
Visualização da contrução do Modelo 2.2 para as posições angulares θ = π/18, Fig. (a), θ = −π/18, Fig. (b), e θ = 0, Fig. (c), onde todas essas posições estão presentes na cena 1 nomeada de Pêndulo físico acoplado. Temos também os gráficos do movimento para a posição angular θ(t), velocidade angular ω(t) e aceleração angular α(t), respectivamente nas cores vermelho, preto e zul, apresentado na Fig. (d) para a cena 2 nomeada de Gráficos para o movimento.

Para a contrução da cena 2, nomeada de Gráficos para o movimento, usamos a função gdots() e a condicional if com a limitação em t ˂ 10 à visualização da animação dos gráficos da posição θ(t), velocidade ω(t) e aceleração α(t), respectivamente nas cores vermelho, preto e zul, como mostrado na cena 2, Fig. 12(d). O código completo utilizado na construção do Modelo 2.2 pode ser visto no Apêndice [codigo2]C.

4.2.3. Construção e visualização do Modelo 2.3

No oscilador acoplado, Modelo 2.3, construímos a animação através de três cenas: cena 1 – visualização da dinâmica do oscilador; cena 2 – apresentação dos gráficos da posição x₁(t), velocidade v₁(t) e aceleração a₁(t) para o objeto 1 (caixa 1); e cena 3 – apresentação dos gráficos da posição x₂(t), velocidade v₂(t) e aceleração a₂(t) para o objeto 2 (caixa 2). A contrução da cena 1 para o oscilador acoplado é baseda na construção do Modelo 2.1, a diferença aqui é que temos duas molas e duas caixas. Para definirmos os objetos usamos as funções box() e helix(), as mesmas já foram discutidas no Modelo 2.1.

Consideramos os seguintes parâmetros em unidades arbitrárias para o modelo: massas m₁ = m₂ = m = 5.0, constantes elásticas das molas k₁ = k₂ = 1.0; condições iniciais x₁(0) = 1.0, v₁(0) = 0.0, x₂(0) = -2.0 e v₂(0) = 0.0; e passo de tempo δt = 0.01. Usamos a função while() em conjunto com as equações de movimento para os dois objetos, são elas: aceleração a₁ = (k₁/m) ∗(x₂−2 ∗x₁) e a₂ = −(k₂/m) ∗ (x₂ −x₁); velocidade v₁ + = a₁ ∗ dt e v₂ + = a₂ ∗ dt; posição x₁ + = v₁ ∗ dt e x₂ + = v₂ ∗ dt. As equações para a₁ e a₂ foram obtidas com a aplicação da Segunda Lei de Newton dadas respectivamente pelas equações (14) e (15).

Ainda na cena 1 temos que a animação gerada para os blocos e as molas foram associadas às equaçãoes de posição. Sendo assim, temos para o movimento do bloco 1 a expressão bloco1.pos.x = x₁ - 3, para o bloco 2 temos a expressão bloco2.pos.x = x₂ + 3 e paras as molas 1 e 2 escrevemos respectivamente as expressões mola1.axis = bloco1.pos.x − parede1.pos.x e mola2.axis = bloco2.pos.x − bloco1.pos.x. Algumas posições específicas para os blocos 1 e 2 são mostradas na Fig. 13(a) para as condições iniciais x₁ = 1 e x₂ = −2, na Fig. 13(b) temos x₁ = 0 e x₂ = 0 que são as posições de equilíbrio do sistema, e na Fig. 13(c) temos x₁ = −1 e x₂ = 2.

Figura 13
Visualização da contrução do Modelo 2.3 para as posições x₁ = 1 e x₂ = −2, Fig. (a); x₁ = 0 e x₂ = 0, Fig. (b); e x₁ = −1 e x₂ = 2, Fig. (c), onde todas essas posições estão presentes na cena 1. Para as cenas 2 e 3, respectivimamente nomeadas de Gráficos para o movimento do bloco 1 e Gráficos para o movimento do bloco 2, temos os gráficos das posições x₁(t) e x₂(t) (curvas em vermelho), velocidades v₁(t) e v₂(t) (curvas em preto) e acelerações a₁(t) e a₂(t) (curvas em zul). A cena 2 é mostrada na Fig. (d) e a cena 3 na Fig. (e).

A construção do Modelo 2.3 é finalizada com as criações das cenas 2 e 3, respectivamente mostradas nas Figs. 13(d) e 13(e). Em ambas as cenas usamos a função condicional if para limitarmos a visualização da animação dos gráficos em t ˂ 40. A cena 2 apresenta os gráficos da posição x₁(t) (curva em vermelho), velocidade v₁(t) (curva em preto) e aceleração a₁(t) (curva em zul) par ao bloco 1, essa cena foi nomeada de Gráficos para o movimento do bloco 1. A cena 3 apresenta os gráficos da posição x₂(t) (curva em vermelho), velocidade v₂(t) (curva em preto) e aceleração a₂(t) (curva em zul) par ao bloco 2, essa cena foi nomeada de Gráficos para o movimento do bloco 2. O código completo utilizado na construção do Modelo 2.3 pode ser visto no Apêndice codigo3[].

5. Conclusão

Neste trabalho, apresentamos os resultados obtidos por meio da comparação entre soluções exatas e soluções aproximadas para problemas de três tipos de osciladores: sistema massa-mola com rotação, pêndulo físico acoplado a uma mola e um sistema de oscilador acoplado. Para o estudo numérico, usamos a linguagem de programação Python e suas bibliotecas NumPy, Pandas e Matplotlib. Mostramos também os resultados obtidos na construção e a visualização da animação dos osciladores, através da biblioteca VPython. Destacamos a importância de incorporar a linguagem de programação Python e suas bibliotecas para solucionarmos numericamente as equações diferenciais de segunda ordem dos três modelos de osciladores, além de visualizar graficamente as soluções analíticas e aproximadas. As soluções aproximadas foram obitidas com a implementação computacional dos métodos numéricos de Euler, de Heun e o método de Runge-Kutta de quarta ordem, onde foi considerando três passos de tempo, δt = t_n+1−t_n, de 0.2 e 0.05 para cada um dos métodos. O comparativo entre soluções exatas e soluções aproximadas nos deu uma compreensão aprofundada e permitiu explorararmos de maneira significativa as soluções econtradas. Obtivemos então que o método de Runge-Kutta de quarta ordem nos forneceu uma melhor aproximação da solução exata nos valores de posição e velocidade dos osciladores. Acreditamos que esse trabalho tem o potencial de ajudar, através de uma abordagem mais interativa, na análise e o entedimento mais completo de problemas de osciladores para serem aplicados nas aulas de física nos cursos de graduação em física, exatas e engenharias como uma alternativa de experimento virtual.

Agradecimentos

Edinardo I. B. Rodrigues e M. L. da S. Dias agradecem ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) e ao Instituto Federal de Educação, Ciência e Tecnologia Baiano (IF Baiano) pelo suporte durante o desenvolvimento deste trabalho.

Apêndice A: Solução Para a Equação de Movimento do Oscilador Acoplado

Considere um sistema livre de atrito onde dois corpos, de massas iguais a m que estão presos por meio de duas molas de constantes elásticas k₁ e k₂, podem mover-se livremente, conforme mostrado na Fig. 3. A força resultante que o corpo 1 sente devido as ações das molas é

(A.1)$$F_1-F_2=-k_1{x}_1+k_2(x_2-x_1)$$

e para o corpo 2 temos que a força resultante é dada por

(A.2)$$F_2=-k_2(x_2-x_1).$$

Aplicando a segunda lei de Newton temos o par de equações diferenciais

(A.3)$$\frac{d^2{x}_1(t)}{d{t}^2}=-\frac{k_1+k_2}{m}{x}_1(t)+\frac{k_2}{m}{x}_2(t)$$

e

(A.4)$$\frac{d^2{x}_2(t)}{d{t}^2}=\frac{k_2}{m}{x}_1(t)-\frac{k_2}{m}{x}_2(t).$$

É razoável assumir que o movimento resultante tem um comportamento oscilatório e como buscamos modos normais, onde todas as partes do sistema acoplado oscilam na mesma frequência, consideramos as seguintes funções de teste para x₁(t) e x₂(t)

(A.5)$$x_1(t)=A{e}^{i\omega t},$$

(A.6)$$x_2(t)=B{e}^{i\omega t},$$

onde A e B são amplitudes a serem determinadas. Substituindo essas funções nas equações (A.1) e (A.4), obtemos que

(A.7)$$\left(-{\omega }^2+\frac{k_1+k_2}{m}\right)A-\frac{k_2}{m}B=0,$$

(A.8)$$-\frac{k_2}{m}A+\left(-{\omega }^2+\frac{k_2}{m}\right)B=0.$$

Podemos reescrever esse sistema na forma matricial, então

(A.9)$$\left[\begin{array}{cc}-{\omega }^2+\frac{k_1+k_2}{m}& -\frac{k_2}{m}\\ -\frac{k_2}{m}& -{\omega }^2+\frac{k_2}{m}\end{array}\right]\left[\begin{array}{c}A\\ B\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

Este sistema tem apenas uma solução não trivial para o determinante da matriz dos coeficientes ser nulo, portanto

(A.10)$$m^2{\omega }^4-m(k_1+2k_2){\omega }^2+k_1{k}_2=0.$$

Resolvendo essa equação obtemos as seguintes soluções

(A.11)$${\omega }_1=\pm {\left[\frac{k_1+2k_2}{2m}-\frac{\sqrt{k_1^2+4k_2^2}}{2m}\right]}^{1/2},$$

(A.12)$${\omega }_2=\pm {\left[\frac{k_1+2k_2}{2m}+\frac{\sqrt{k_1^2+4k_2^2}}{2m}\right]}^{1/2}.$$

A relação entre os coeficientes A e B pode ser obitida da equação (A.7), então

(A.13)$$B=\frac{-m{\omega }^2+k_1+k_2}{k_2}A.$$

Combinando essa relação com ω₁, equação (A.11), temos que

(A.14)$$B_1=\frac{A_1}{2k_2}\left(k_1+\sqrt{k_1^2+4k_2^2}\right),$$

e combinando (A.13) com ω₁, equação (A.12), obtemos

(A.15)$$B_2=\frac{A_2}{2k_2}\left(k_1-\sqrt{k_1^2+4k_2^2}\right).$$

Uma vez obitidas as frequências ω e ω₂ e os coeficientes B₁ e B₂, podemos escrever que a solução mais geral, x₁(t) e x₂(t), para as equações de movimento, equações (A.3) e (A.4), podem ser escritas nas formas abaixo:

(A.16)$$x_1(t)=A_1\text{cos}({\omega }_{1}t)+A_2\text{cos}({\omega }_{2}t)$$

e

(A.17)$$x_2(t)=B_1\text{cos}({\omega }_{1}t)+B_2\text{cos}({\omega }_{2}t).$$

Apêndice B: Código Para a Construção da Visualização da Animação do Modelo 2.1

Apêndice C: Código Para a Construção da Visualização da Animação do Modelo 2.2

Apêndice D: Código Para a Construção da Visualização da Animação do Modelo 2.3

Referencias

• ^[1] F.A. Moura, P.H.P. Mandarino e S.C.P. Silva, Rev. Bras. Ens. Fis. 42, e20200044 (2020).
• ^[2] V.A. Oliveira e A.C. da Silva, Ens. Pesqui. Educ. Ciênc. 22, e11969 (2020).
• ^[3] M.A. Moreira, Rev. Bras. Ens. Fis. 43, e20200451 (2021).
• ^[4] R.C. Barbosa, J.F. Junior e P.L. Junior, Rev. Bras. Ens. Fis. 45, e20230210 (2023).
• ^[5] N.T.O. Reis, N.M.D. Garcia, P.N. Souza e P.S. Baldessar, Rev. Bras. Ens. Fis. 30, 1401 (2008).
• ^[6] C.H.J. Davies, Computers & Education 39, 271 (2002).
• ^[7] M.A. Cavalcante e C.R.C. Tavolaro, Rev. Bras. Ens. Fis. 22, 421 (2000).
• ^[8] G. Camiletti e L.A. Ferracioli, Rev. Bras. Ens. Fis. 24, 110 (2002).
• ^[9] G.R. Neves, J.C. Silva e V.D. Teodoro, arXiv:1006.4662v1 (2010).
• ^[10] D.M. Jaime e A.A. Leonel, Rev. Bras. Ens. Fis. 46, 20230309 (2024).
• ^[11] J.A. Macêdo, A.G. Dickman e I.S.F. Andrade, Cad. Bras. de Ens. de Fís. 29, 562 (2012).
• ^[12] L.R. Petzold, L.O. Jay e J. Yen, Ada Numerica 6, 437 (1997).
• ^[13] J. Bartos, em: In WDS’05 Proceedings of Contributed Papers, Part III (Praga, 2005).
• ^[14]D. Schott., Glob. Jour. of Eng. Educ. 12, 17 (2010).
• ^[15] L.J. Kwari, J. Sunday, J.N. Ndam e A.A. James, Science Forum (Jour. of Pure and Appl. Sc.) 21, 503 (2021).
• ^[16] N.M. Janah, F. Falah, Ratnawati, Ishafit e W.S.B. Dwandaru, Ind. Rev. of Phys. 4, 46 (2021).
• ^[17] E.I.B. Rodrigues, Rev. Velho Chico 2, 15 (2022).
• ^[18] S.V. Kontomaris, G. Chliveros e A. Malamou, Mult. Scient. Jour. 6, 140 (2023).
• ^[19] C.P.D. Conceição e T.D. Admiral, Braz. Jour. of Science 1, 64 (2022).
• ^[20] L.E. Borges, Python para Desenvolvedores (Novatec, Rio de Janeiro, 2014), 3 ed.
• ^[21] L. Euler, Institutionum calculi integralis volumen primum (Académie Impériale des Sciences, Saint-Petersbourg, 1768).
• ^[22] A.L. Alves,S.S. Bento e C.H. Marchi, Rev. Bras. Ens. Fis. 45, e20230229 (2023).
• ^[23] K. Heun e Z. Math. Phys. 45, 23 (1900).
• ^[24] C. Runge, Mathematische Annalen 46, 167 (1895).
• ^[25] Kutta, Zeitschrift für Mathematik und Physik 46, 435 (1901).
• ^[26] R.L. Sterza e A.C. Brandi, Rev. Eletr. Paul. de Matem. 7, 12 (2016).
• ^[27] K.A. Koroche, Int. Jour. of Syst. Sc. and Appl. Math. 6, 1 (2021).
• ^[28] P.F.T. Dorneles, I.S. Araujo e E.A. Veit, Rev. Bras. Ens. Fis. 30, 3 (2008).
• ^[29] R.L. Viana, Rev. Bras. Ens. Fís. 33, 3 (2011).
• ^[30] A. Kandus, F.W. Gutmann e C.M.C. de Castilho, Rev. Bras. Ens. Fis. 28, 4 (2006).
• ^[31] F.F. Pedrozo e T.C. de Freitas, Rev. Bras. Ens. Fis. 44, (2022).
• ^[32] J. Jani, WSEAS Transactions on Advances in Engineering Education 20, 1 (2023).
• ^[33] G. Akpojotor, L. Ehwerhemuepha, M. Echenim e F. Akpojotor, African Journal of Physics 3, 94 (2010).
• ^[34] E. Tufino, S. Oss e L. Gratton, J. Phys.: Conf. Ser. 2490, 012003 (2023).
• ^[35] P.D. Bhoyar, Modern trends in multidisciplinary subjects (Redshie Publication, Lunawada, 2023), v. 3.
• ^[36] K. Shah, TechRxiv:10.36227/techrxiv.22581253.v1 (2023).
• ^[37] S.R. de Oliveira, Rev. Bras. Ens. Fis. 46, (2024).
• ^[38] NumPy user guide, disponível em:https://numpy.org, acessado em: 10/06/2024.
  » https://numpy.org
• ^[39] User Guide, disponível em:https://pandas.pydata.org, acessado em: 10/06/2024.
  » https://pandas.pydata.org
• ^[40] Using Matplotlib, disponível em:https://matplotlib.org/, acessado em: 10/06/2024.
  » https://matplotlib.org/
• ^[41] Python, disponível em:https://www.python.org, acessado em: 15/06/2024.
  » https://www.python.org
• ^[42] VPython 3D Programming for Ordinary Mortals, disponível em:https://www.vpython.org/index.html, acessado em: 15/06/2024.
  » https://www.vpython.org/index.html
• ^[43] F.L.Q. Vidal e C.A.R.S. Fontoura, Open Science Research 1, 1108 (2022).
• ^[44] C.S. Cruz, L.Q. Galvão, S. Rosa e W.S. Santana, Cad. Bras. de Ens. De Fís. 9, 204 (2022).
• ^[45] VPython documentation, disponível em:https://www.glowscript.org/docs/VPythonDocs/index.html, acessado em: 16/06/2024.
  » https://www.glowscript.org/docs/VPythonDocs/index.html

Editado por

• Editor-Chefe:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Datas de Publicação

• Publicação nesta coleção
  31 Mar 2025
• Data do Fascículo
  2025

Histórico

• Recebido
  25 Nov 2024
• Revisado
  12 Fev 2025
• Aceito
  18 Fev 2025