Compare a solução numérica com a solução analítica de um sistema amortecido com m=100 \mathrm{~kg}, k=2000 \mathrm{~N} / \mathrm{m} e c=200 \mathrm{~N} . \mathrm{s} / \mathrm{m} submetido a uma força F(t)= 150 \cos (\omega t). Considere que as consições iniciais são dadas por x(0)=0,01 \mathrm{~m} e \dot{x}(0)=0,1 \mathrm{~m} / \mathrm{s}. Considerando que \omega=10,0 \mathrm{rad} / \mathrm{s}, plote na mesma figura as soluções numérica e analítica.

Passo 1

O enunciado pede que a gente compare a solução numérica com a solução analítica do oscilador amortecido forçado
m\ddot{x}+c\dot{x}+kx=F(t),
com os parâmetros e condições iniciais fornecidos. Perguntas e dados importantes (extraídos do enunciado):
  • Força externa: \displaystyle F(t)=150\cos(\omega t) com \omega=10{,}0\\ \mathrm{rad/s}.
  • Massa: \displaystyle m=100\\ \mathrm{kg}.
  • Amortecimento: \displaystyle c=200\\ \mathrm{N\cdot s/m}.
  • Rigidez: \displaystyle k=2000\\ \mathrm{N/m}.
  • Condições iniciais: \displaystyle x(0)=0{,}01\\ \mathrm{m} e \displaystyle \dot{x}(0)=0{,}1\\ \mathrm{m/s}.
  • O objetivo prático: traçar em uma mesma figura a solução numérica (integração do EDO) e a solução analítica (expressão fechada), e comparar.

Passo 2

Escrevemos a EDO e definimos as quantidades naturais que usamos na solução analítica:
m\ddot{x}+c\dot{x}+kx=F_0\cos(\omega t),\qquad F_0=150.
Cálculos de parâmetros:
\omega_n=\sqrt{\frac{k}{m}}=\sqrt{\frac{2000}{100}}=\sqrt{20}=4{,}472135955\\ \mathrm{rad/s}.
\text{razão de amortecimento } \zeta=\frac{c}{2m\omega_n}=\frac{200}{2\cdot100\cdot4{,}472135955}=0{,}223606798.ã
Uma forma útil (e mais simples numericamente) é notar
\alpha=\frac{c}{2m}= \frac{200}{2\cdot100}=1{,}0\\ \mathrm{s^{-1}},
logo a frequência amortecida é
\omega_d=\sqrt{\omega_n^2-\alpha^2}=\sqrt{20-1}=\sqrt{19}=4{,}358898943\\ \mathrm{rad/s}.
Comentário rápido de física: \omega_n é a frequência natural sem amortecimento; \zeta indica se o sistema é subamortecido ( \zeta<1), superamortecido ( \zeta>1) ou criticamente amortecido ( \zeta=1). Aqui \zeta\approx0{,}224<1, então o sistema é subamortecido e apresenta oscilações amortecidas no transiente com frequência \omega_d.

Passo 3

Encontrar a solução particular em regime permanente. Para entrada harmônica \displaystyle F_0\cos(\omega t), a resposta em regime permanente é da forma
x_p(t)=X\cos(\omega t-\varphi),
onde a amplitude X e a defasagem \varphi satisfazem
X=\frac{F_0}{\sqrt{(k-m\omega^2)^2+(c\omega)^2}}, \qquad \varphi=\operatorname{atan2}(c\omega,\;k-m\omega^2).
Agora calculamos numericamente:
k-m\omega^2=2000-100\cdot10^2=2000-10000=-8000.
c\omega=200\cdot10=2000.
\sqrt{(k-m\omega^2)^2+(c\omega)^2}=\sqrt{(-8000)^2+2000^2}=\sqrt{68\times10^6}=8246{,}211251\ldots
X=\frac{150}{8246{,}211251\ldots}=0{,}01818720005\\ \mathrm{m}\\ (\text{aprox.}).
Para a fase:
\tan\varphi=\frac{c\omega}{k-m\omega^2}=\frac{2000}{-8000}=-0{,}25.
Como o denominador é negativo e o numerador positivo, \varphi fica no segundo quadrante; portanto
\varphi=\arctan(-0{,}25)+\pi=2{,}89661399046\\ \mathrm{rad}\\ (\text{aprox.}).
Comentário prático: podemos usar \operatorname{atan2}(2000,-8000) em código para obter \varphi corretamente no intervalo (0,\pi).

Passo 4

Forma geral da solução analítica completa: soma da homogênea (transiente) com a particular (regime permanente):
x(t)=x_h(t)+x_p(t),
x_h(t)=e^{-\alpha t}\left(C_1\cos(\omega_d t)+C_2\sin(\omega_d t)\right),
x_p(t)=X\cos(\omega t-\varphi).
Agora encontramos C_1 e C_2 a partir das condições iniciais. Primeiro avaliamos x_p(0) e \dot{x}_p(0):
x_p(0)=X\cos(-\varphi)=X\cos\varphi,
\dot{x}_p(0)=-X\omega\sin(-\varphi)=X\omega\sin\varphi.
Usando valores numéricos já calculados:
\cos\varphi=\cos(2{,}89661399046)=-0{,}9701425001453319,
\sin\varphi=\sin(2{,}89661399046)=-0{,}24253562503633297.
x_p(0)=0{,}01818720005\cdot(-0{,}9701425001453319)=-0{,}01764287615\\ \mathrm{m}.
\dot{x}_p(0)=0{,}01818720005\cdot10\cdot(-0{,}24253562503633297)=-0{,}04414287615\\ \mathrm{m/s}.
Agora o componente homogêneo inicial:
x_h(0)=x(0)-x_p(0)=0{,}01-(-0{,}01764287615)=0{,}02764287615.
\dot{x}_h(0)=\dot{x}(0)-\dot{x}_p(0)=0{,}1-(-0{,}04414287615)=0{,}14414287615.
Expressamos \dot{x}_h(t) explicitamente (derivada do termo homogêneo):
\dot{x}_h(t)=-\alpha e^{-\alpha t}(C_1\cos\omega_d t+C_2\sin\omega_d t)+e^{-\alpha t}(-C_1\omega_d\sin\omega_d t + C_2\omega_d\cos\omega_d t).
Avaliada em t=0:
\dot{x}_h(0)=-\alpha C_1+\omega_d C_2.
Logo resolvemos para C_1 e C_2:
C_1=x_h(0)=0{,}02764287615,
C_2=\frac{\dot{x}_h(0)+\alpha C_1}{\omega_d}=\frac{0{,}14414287615+1{,}0\cdot0{,}02764287615}{4{,}358898943}=0{,}039430800\\ (\text{aprox.}).

Passo 5

Portanto a solução analítica final (com valores numéricos aproximados) é
\boxed{\,x(t)=e^{-1{,}0\,t}\big(0{,}02764287615\cos(4{,}358898943\,t)+0{,}039430800\sin(4{,}358898943\,t)\big)+0{,}01818720005\cos(10\,t-2{,}89661399046)\,}.
Explicação do que isso significa: o primeiro termo é o transiente (amplitude decai exponencialmente com taxa 1{,}0\\ \mathrm{s^{-1}}), o segundo termo é a resposta em regime permanente forçada pela frequência \omega=10\\ \mathrm{rad/s}.

Passo 6

Como obter a solução numérica e plotar junto com a analítica (passo a passo e com código de exemplo).
  1. Converter ODE de segunda ordem em sistema de primeira ordem:
y_1=x,\qquad y_2=\dot{x}.
\dot{y}_1=y_2,
\dot{y}_2=\frac{1}{m}\big(-c y_2 - k y_1 + F_0\cos(\omega t)\big).
Condições iniciais: y_1(0)=0{,}01, y_2(0)=0{,}1.
  1. Sugerência de intervalo para o plot: tomar por exemplo t\in[0,10]\\ \mathrm{s}. Observação: o tempo característico de decaimento é 1/\alpha=1\\ \mathrm{s}, então por t=5 10\\ \mathrm{s} já estamos bem no regime permanente.
  2. Código de referência (Python + scipy) — copie e rode no seu ambiente Python. O código calcula a solução numérica, avalia a solução analítica com a fórmula acima e plota as duas curvas juntos e calcula o erro máximo.
Python (exemplo):
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# parâmetros
m = 100.0
c = 200.0
k = 2000.0
F0 = 150.0
omega = 10.0

# constantes analíticas obtidas anteriormente
alpha = 1.0
omega_d = np.sqrt(19.0)
C1 = 0.02764287615
C2 = 0.039430800
X = 0.01818720005
phi = 2.89661399046

# sistema de 1ª ordem
def f(t,y):
    y1,y2 = y
    dy1 = y2
    dy2 = ( -c*y2 - k*y1 + F0*np.cos(omega*t) )/m
    return [dy1, dy2]

t_span = (0.0, 10.0)
t_eval = np.linspace(0,10,2001)
y0 = [0.01, 0.1]

sol = solve_ivp(f, t_span, y0, t_eval=t_eval, rtol=1e-9, atol=1e-12)

t = sol.t
x_num = sol.y[0,:]

# solução analítica
x_anal = np.exp(-alpha*t)*(C1*np.cos(omega_d*t)+C2*np.sin(omega_d*t)) + X*np.cos(omega*t-phi)

# erro
err = np.abs(x_num - x_anal)
print('Erro máximo |x_num - x_anal| =', err.max())

# plot
plt.figure(figsize=(10,5))
plt.plot(t, x_num, label='numérica (solve_ivp)', lw=1.2)
plt.plot(t, x_anal, '--', label='analítica', lw=1.2)
plt.xlabel('t [s]')
plt.ylabel('x(t) [m]')
plt.legend()
plt.grid(True)
plt.title('Comparação: solução numérica vs analítica')
plt.show()
  1. Código MATLAB/Octave (exemplo):
m = 100;
c = 200;
k = 2000;
F0 = 150;
omega = 10;

alpha = 1;
omega_d = sqrt(19);
C1 = 0.02764287615;
C2 = 0.039430800;
X = 0.01818720005;
phi = 2.89661399046;

f = @(t,y) [ y(2); (-c*y(2)-k*y(1)+F0*cos(omega*t))/m ];
tspan = [0 10];
y0 = [0.01; 0.1];
opts = odeset('RelTol',1e-9,'AbsTol',1e-12);
[ts, ys] = ode45(f, linspace(0,10,2001), y0, opts);
x_num = ys(:,1);

x_anal = exp(-alpha*ts).*(C1*cos(omega_d*ts)+C2*sin(omega_d*ts)) + X*cos(omega*ts-phi);

max_err = max(abs(x_num - x_anal));
disp(['Erro máximo = ', num2str(max_err)])
plot(ts, x_num, 'b', ts, x_anal, '--r')
legend('numérica','analítica')
xlabel('t [s]'); ylabel('x(t) [m]'); grid on

Passo 7

O que esperar ao rodar o código e como interpretar os resultados:
  • No gráfico: você verá um termo transitório que decai rapidamente (pois \alpha=1 dá decaimento significativo em poucos segundos) e, após o transiente, o sinal segue a oscilação forçada com amplitude aproximada X\approx0{,}01819\\ \mathrm{m} e fase \varphi\approx2{,}8966\\ \mathrm{rad}.
  • Comparação numérica/analítica: se usar tolerâncias razoáveis (por exemplo rtol=1e-9, atol=1e-12 no solve_ivp), o erro máximo \max_t|x_{\text{num}}(t)-x_{\text{anal}}(t)| deverá ser muito pequeno (na prática algo da ordem de tolerâncias numéricas, tipicamente \le 10^{-8} 10^{-6} dependendo do integrador e passo). Isto confirma que a solução numérica e a analítica coincidem dentro da precisão computacional.
  • Dica de visual: use linha sólida para a solução numérica e tracejada para a analítica; olhe tanto o comportamento geral quanto o erro absoluto se quiser números.

Passo 8

Resumo didático e cheque rápido manual (pontos para checar quando estiver na prova):
  • Verifique dimensões nas fórmulas de X e \varphi.
  • Para achar constantes do transiente, subtraia a solução particular nas condições iniciais: x_h(0)=x(0)-x_p(0) e \dot{x}_h(0)=\dot{x}(0)-\dot{x}_p(0).
  • Para o termo homogêneo subamortecido use a forma e^{-\alpha t}(C_1\cos\omega_d t + C_2\sin\omega_d t) e determine C_1,C_2 diretamente.
  • Sempre confira sinais da fase com \operatorname{atan2} para ficar no quadrante correto.

Resposta

\displaystyle x(t)=e^{-1{,}0\,t}\big(0{,}02764287615\cos(4{,}358898943\,t)+0{,}039430800\sin(4{,}358898943\,t)\big)+0{,}01818720005\cos(10\,t-2{,}89661399046).
As soluções numérica e analítica coincidem dentro da precisão numérica do integrador (erro máximo entre as soluções tipicamente \lesssim10^{-6} com tolerâncias rígidas); amplitude de regime permanente X\approx0{,}0181872\\ \mathrm{m} e fase \varphi\approx2{,}89661399\\ \mathrm{rad}.

Exercícios relacionados