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:
  • Massa: m=10\\ \mathrm{kg}.
  • Rigidez da mola: k=400\\ \mathrm{N/m}.
  • Força externa: F(t)=400\cos(\omega t). Portanto F_0=400.
  • 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
\omega_n^2 = \frac{k}{m}.
Substituímos os valores dados:
\omega_n^2 = \frac{400}{10} = 40,\qquad \omega_n = \sqrt{40} = 2\sqrt{10}\approx 6{,}3245553203.
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):
  1. 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.
  2. 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)
(-\omega^2 + \omega_n^2)C\cos(\omega t) = \frac{F_0}{m}\cos(\omega t).
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.
  1. A solução geral é
x(t) = A\cos(\omega_n t) + B\sin(\omega_n t) + C\cos(\omega t).
Impondo as condições iniciais x(0)=0 e \dot{x}(0)=10:
  • x(0)=0 A + C = 0 \Rightarrow A=-C.
  • Derivando,
\dot{x}(t) = -A\omega_n\sin(\omega_n t) + B\omega_n\cos(\omega_n t) - C\omega\sin(\omega t).
Em t=0: \dot{x}(0)=B\omega_n = 10 \Rightarrow B = \dfrac{10}{\omega_n}.
  1. Substituindo A=-C e B na solução geral, obtemos uma forma bem prática e compacta:
x_{\text{an}}(t) = C\bigl(\cos(\omega t) - \cos(\omega_n t)\bigr) + \frac{10}{\omega_n}\sin(\omega_n t),
com
C = \frac{F_0}{k - m\omega^2}.
Essa expressão é a nossa solução analítica explícita (válida para \omega\neq\omega_n).

Passo 4

Aplicando a fórmula para os valores numéricos (calculei os coeficientes para usar nos plots também):
  • Valor comum:
\frac{10}{\omega_n} = \frac{10}{\sqrt{40}} = \frac{\sqrt{10}}{2}\approx 1{,}5811388301.
  • Para \omega=20{,}0:
k - m\omega^2 = 400 - 10\cdot 20^2 = 400 - 4000 = -3600,
logo
C = \frac{400}{-3600} = -\frac{1}{9} \approx -0{,}1111111111.
  • Para \omega=20{,}1:
k - m\omega^2 = 400 - 10\cdot (20{,}1)^2 = 400 - 4040{,}1 = -3640{,}1,
logo
C \approx \frac{400}{-3640{,}1} \approx -0{,}109890\\ (\text{aprox.}).
  • Para \omega=35{,}0:
k - m\omega^2 = 400 - 10\cdot 35^2 = 400 - 12250 = -11850,
logo
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):
  1. 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.
  1. 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).
  2. 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
\text{err}_{\infty} = \max_{t} |x_{\text{num}}(t) - x_{\text{an}}(t)|
e o erro RMS
\text{err}_{\mathrm{RMS}} = \sqrt{\frac{1}{N}\sum_{i=1}^N (x_{\text{num}}(t_i)-x_{\text{an}}(t_i))^2 }.]

Passo 6

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)
def forced_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ção
def x_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étricas
    print(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.

Resposta

x_{\text{an}}(t)\;=\;C\bigl(\cos(\omega t)-\cos(\sqrt{40}\,t)\bigr)+\frac{10}{\sqrt{40}}\sin(\sqrt{40}\,t),
onde C=\dfrac{400}{400-10\omega^2}=\dfrac{400}{k-m\omega^2}.
Para os três casos:
\text{(i) }\omega=20{,}0:\qquad C=-\frac{1}{9},\quad x_{\text{an}}(t)=-\frac{1}{9}\bigl(\cos(20t)-\cos(\sqrt{40}\,t)\bigr)+\frac{\sqrt{10}}{2}\sin(\sqrt{40}\,t).
\text{(ii) }\omega=20{,}1:\qquad C\approx -0{,}109890,\quad x_{\text{an}}(t)\approx -0{,}109890\bigl(\cos(20{,}1\,t)-\cos(\sqrt{40}\,t)\bigr)+\frac{\sqrt{10}}{2}\sin(\sqrt{40}\,t).
\text{(iii) }\omega=35{,}0:\qquad C=-\frac{8}{237}\approx -0{,}033755,\quad x_{\text{an}}(t)\approx -0{,}033755\bigl(\cos(35t)-\cos(\sqrt{40}\,t)\bigr)+\frac{\sqrt{10}}{2}\sin(\sqrt{40}\,t).
(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.)

Exercícios relacionados