Compare a solução numérica com a solução analítica de um sistema amortecido com
e
submetido a uma força
. Considere que as consições iniciais são dadas por
e
. Considerando que
, 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
com os parâmetros e condições iniciais fornecidos. Perguntas e dados importantes (extraídos do enunciado):
- Força externa:
com . - Massa:
. - Amortecimento:
. - Rigidez:
. - Condições iniciais:
e . - 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:
Cálculos de parâmetros:
Uma forma útil (e mais simples numericamente) é notar
logo a frequência amortecida é
Comentário rápido de física:
é a frequência natural sem amortecimento;
indica se o sistema é subamortecido (
), superamortecido (
) ou criticamente amortecido (
). Aqui
, então o sistema é subamortecido e apresenta oscilações amortecidas no transiente com frequência
.
Passo 3
Encontrar a solução particular em regime permanente. Para entrada harmônica
, a resposta em regime permanente é da forma
onde a amplitude
e a defasagem
satisfazem
Agora calculamos numericamente:
Para a fase:
Como o denominador é negativo e o numerador positivo,
fica no segundo quadrante; portanto
Comentário prático: podemos usar
em código para obter
corretamente no intervalo
.
Passo 4
Forma geral da solução analítica completa: soma da homogênea (transiente) com a particular (regime permanente):
Agora encontramos
e
a partir das condições iniciais. Primeiro avaliamos
e
:
Usando valores numéricos já calculados:
Agora o componente homogêneo inicial:
Expressamos
explicitamente (derivada do termo homogêneo):
Avaliada em
:
Logo resolvemos para
e
:
Passo 5
Portanto a solução analítica final (com valores numéricos aproximados) é
Explicação do que isso significa: o primeiro termo é o transiente (amplitude decai exponencialmente com taxa
), o segundo termo é a resposta em regime permanente forçada pela frequência
.
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).
- Converter ODE de segunda ordem em sistema de primeira ordem:
Condições iniciais:
,
.
-
Sugerência de intervalo para o plot: tomar por exemplo
. Observação: o tempo característico de decaimento é , então por – já estamos bem no regime permanente. -
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()
- 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
dá decaimento significativo em poucos segundos) e, após o transiente, o sinal segue a oscilação forçada com amplitude aproximada e fase . - 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
deverá ser muito pequeno (na prática algo da ordem de tolerâncias numéricas, tipicamente – 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
e . - Para achar constantes do transiente, subtraia a solução particular nas condições iniciais:
e . - Para o termo homogêneo subamortecido use a forma
e determine diretamente. - Sempre confira sinais da fase com
para ficar no quadrante correto.
Resposta
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