Compare a solução numérica com a solução analítica de um sistema massa-mola com m=10 \mathrm{~kg} e k=400 \mathrm{~N} / \mathrm{m} submetido a uma força F(t)=400 \cos (\omega t). Considere que as consições iniciais são dadas por x(0)=0 e \dot{x}(0)=10 \mathrm{~m} / \mathrm{s}. Considerando que \omega=20,0 \mathrm{rad} / \mathrm{s}, plote na mesma figura as soluções numérica e analítica. Refaça o exercício \operatorname{com} \omega=20,1 \mathrm{rad} / \mathrm{s} e \omega=35,0 \mathrm{rad} / \mathrm{s}.
Passo 1
O enunciado pede que a gente compare a solução numérica e a solução analítica do seguinte sistema massa-mola forçado, com condições iniciais dadas, para três valores de \omega. As informações importantes que temos são:
Condições iniciais: x(0)=0 e \dot{x}(0)=10\\ \mathrm{m/s}.
Três casos para comparar:
\omega=20{,}0\\ \mathrm{rad/s}.
\omega=20{,}1\\ \mathrm{rad/s}.
\omega=35{,}0\\ \mathrm{rad/s}.
Objetivo: obter a solução analítica x_{\text{an}}(t), obter solução numérica x_{\text{num}}(t) (por integração numérica do ODE), e plotar as duas na mesma figura para cada valor de \omega. Além disso é desejável quantificar brevemente a diferença (erro) entre as soluções.
Passo 2
Escrevendo o problema matematicamente: a equação do movimento é
m\ddot{x} + k x = F_0\cos(\omega t).
Dividindo por m,
\ddot{x} + \omega_n^2 x = \frac{F_0}{m}\cos(\omega t),
onde definimos a frequência natural não amortecida
Vou usar \omega_n=\sqrt{40} nas fórmulas exatas quando apropriado.
Passo 3
Método para solução analítica (passo-a-passo didático):
Ideia: a equação é linear com coeficientes constantes e termo forçante harmónico. A solução geral é soma da solução homogênea (complementar) e de uma solução particular.
Solução homogênea da equação
\ddot{x} + \omega_n^2 x = 0
é
x_h(t) = A\cos(\omega_n t) + B\sin(\omega_n t),
onde A e B são constantes determinadas pelas condições iniciais.
3) Procuramos uma solução particular do tipo harmônico: suponha
x_p(t) = C\cos(\omega t).
Substituindo em \ddot{x} + \omega_n^2 x = \dfrac{F_0}{m}\cos(\omega t) dá
Como \cos(\omega t) não é a função zero identicamente, obtemos (para \omega\neq\omega_n)
C = \frac{F_0/m}{\omega_n^2 - \omega^2} = \frac{F_0}{k - m\omega^2}.
Observação importante: se \omega=\omega_n haveria ressonância e a solução particular seria do tipo t\sin(\omega t). Em nossos três casos \omega não é igual a \omega_n (pois \omega_n\approx 6{,}3246), então a fórmula acima vale.
C = \frac{400}{-11850} = -\frac{8}{237}\approx -0{,}033755.
Eu mantenho as formas exatas nas fórmulas e também os aproximados para implementação numérica.
Passo 5
Como obter a solução numérica (passo-a-passo didático):
Transformar a EDO de segunda ordem em um sistema de primeira ordem. Defina \,y_1=x\) e \(\,y_2=\dot{x}\). Então
\[\begin{cases} \dot{y}_1 = y_2,\\[4pt] \dot{y}_2 = -\omega_n^2 y_1 + \dfrac{F_0}{m}\cos(\omega t). \end{cases}
Condições iniciais: y_1(0)=0\], \(y_2(0)=10.
Escolha de método numérico: um método Runge–Kutta de ordem 4/5 (por exemplo o solver padrão ‘RK45’ do solve_ivp do SciPy) funciona muito bem para este problema não rígido. Alternativamente, para mais precisão ou passo fixo podemos usar um integrador de passo pequeno (RK4 explícito). O importante é escolher um passo de amostragem (para comparar/plotar) suficientemente pequeno para captar a oscilação de maior frequência. A maior frequência angular entre os três casos é \omega=35\\ \mathrm{rad/s}, que corresponde a frequência f=\omega/(2\pi)\approx 5{,}57\\ \mathrm{Hz} e período T\approx 0{,}1796\\ \mathrm{s}. Uma regra prática é tomar pelo menos 100 pontos por período para garantir boa resolução visual e precisão, então \Delta t\lesssim T/100\approx 0{,}0018\\ \mathrm{s}. Vou sugerir \Delta t=0{,}001\\ \mathrm{s} e intervalo de tempo t\in[0,\,5]\\ \mathrm{s} (ou maior se quiser ver mais ciclos).
Implementação numérica: integramos o sistema para obter y_1(t) e comparamos com x_{\text{an}}(t) calculada diretamente pela fórmula do Passo 3. Para medir a discrepância usamos, por exemplo, o erro máximo absoluto
Código Python completo (copie e execute no seu ambiente; ele gera os três gráficos com numérico vs analítico e imprime erros). Comentários no código explicam cada passo.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Parâmetros
m = 10.0
k = 400.0
F0 = 400.0
omega_n = np.sqrt(k/m) # sqrt(40)
v0 = 10.0
x0 = 0.0# Função força (para uso numérico)defforced_system(t, y, omega):
# y[0]=x, y[1]=v
dydt0 = y[1]
dydt1 = - (k/m) * y[0] + (F0/m) * np.cos(omega * t)
return [dydt0, dydt1]
# Solução analítica conforme derivaçãodefx_analytic(t, omega):
C = F0 / (k - m * omega**2)
return C*(np.cos(omega*t) - np.cos(omega_n*t)) + (v0/omega_n)*np.sin(omega_n*t)
# Parâmetros temporais para integração e plot
t0 = 0.0
t_end = 5.0
dt = 0.001
t_eval = np.arange(t0, t_end+dt, dt)
# Valores de omega a estudar
omegas = [20.0, 20.1, 35.0]
for omega in omegas:
# Integração numérica do sistema
sol = solve_ivp(fun=lambda t,y: forced_system(t,y,omega),
t_span=(t0, t_end), y0=[x0, v0], t_eval=t_eval, method='RK45', rtol=1e-9, atol=1e-12)
x_num = sol.y[0]
# Solução analítica no mesmo grid de tempos
x_an = x_analytic(t_eval, omega)
# Erros
err_inf = np.max(np.abs(x_num - x_an))
err_rms = np.sqrt(np.mean((x_num - x_an)**2))
# Imprime métricasprint(f"omega = {omega:.3f} rad/s: max abs error = {err_inf:.3e}, RMS error = {err_rms:.3e}")
# Plot comparativo
plt.figure(figsize=(8,4.5))
plt.plot(t_eval, x_num, label='numérica (RK45)', linewidth=1)
plt.plot(t_eval, x_an, '--', label='analítica', linewidth=1)
plt.xlabel('t (s)')
plt.ylabel('x (m)')
plt.title(f'Solução numérica x_an(x) vs analítica para omega = {omega:.3f} rad/s')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
Explicações sobre escolhas do código:
Usei solve_ivp com tolerâncias pequenas (rtol=1e-9, atol=1e-12) para garantir que o integrador entregue uma solução numérica precisa, assim a discrepância com a solução analítica deve ser limitada à erros numéricos muito pequenos (máquina/integrador).
t_eval é o vetor de tempos em que pedimos o valor da solução; usar passo fixo \Delta t=0{,}001 garante boa resolução temporal para \omega=35\\ \mathrm{rad/s}.
O plot mostra a solução numérica sólida e a analítica tracejada; os erros calculados e impressos servem para comparar numericamente.
Passo 7
Interpretação prática e dicas para a prova:
A derivação analítica é curta e útil: lembrando que a solução particular para F_0\cos(\omega t) assume a forma C\cos(\omega t) com C dado por \dfrac{F_0}{k - m\omega^2} (exceto no caso de ressonância). Em prova, escreva a dedução substituindo x_p=C\cos(\omega t) e resolvendo para C.
Para as condições iniciais, lembro: avalie em t=0 para obter relações entre constantes; derive e avalie em t=0 para obter a outra constante. Isso é padrão e rápido.
Para o método numérico em prova (dependendo do nível exigido), basta mostrar a conversão para sistema de 1ª ordem:\dot{y}_1=y_2,\qquad \dot{y}_2=-\omega_n^2 y_1 + \frac{F_0}{m}\cos(\omega t),e dizer que se usaria RK4 (ou solver automático) com passo \Delta t adequado.
Escolha de \Delta t: sempre garanta várias amostras por período da maior frequência presente no problema. Uma regra prática: pelo menos 50–200 pontos por período.
(Use o código fornecido no Passo 6 para calcular numericamente x_{\text{num}}(t), gerar os gráficos com ambas as soluções sobrepostas e imprimir os erros numéricos.)