<jupyter_start><jupyter_text>Vibrações com excitação harmônica - 2 Graus de Liberdade _Prof. <NAME>_ Considerando o arranjo de massas, molas e amortecedores abaixo:![image.png](attachment:564420c8-601a-4b0a-bc4f-26cf4584d8ba.png)Obtemos o modelo a partir dos diagramas de corpo livre:![image.png](attachment:a71be58a-0407-41ce-92ad-d38b4db5a69c.png)![image.png](attachment:7b2a4552-b2ea-4d4e-8ea2-aff01d9df8ce.png) Para a solução numérica do problema devemos transformar o sistema de duas equações de 2ª Ordem em um sistema de 4 equações de 1ª Ordem.Para isso, utilizamos duas variáveis auxiliares fazendo $y_1 = \dot{x}_1$ e  $y_2 = \dot{x}_2$, e obtemos o seguinte sistema de equações:$\dot{x}_1 = y_1$$\dot{y}_1 = [f_1 -(c_1+c_2)y_1 +c_2*y_2- (k_1+k_2) x_1 + k_2 x_2]/m_1 $$\dot{x}_2 = y_2$$\dot{y}_2 = [f_2 + c_2*y_1 - (c_2+c_3)y_2 + k_2 x_1 - (k_2+k_3) x_2 ]/m_2 $ EXEMPLO: Encontre resposta de vibração para o sistema onde $k_1 = 30$, $k_2 = 5$, $k_3 = 0$, $m_1 = 10 $, $m_2 = 1$, e $c_1 = c_2 = c_3 = 0$ para as condições iniciais $x_1(0) = 1$, $y_1(0) = x_2(0) = y_2(0) = 0$. Para o programa, vamos precisar do módulo  _scipy_ para a integração e o _pylab_ para plotar os resultados.Inicialmente, vamos definir, como uma função, o sistema de equações que será chamado posteriormente:<jupyter_code>def sistema(w, t, p):
    import numpy as nps
    """
    Define o sistema de equações diferenciais para o sistema.

    Argumentos:
        w :  vetor das variáveis de estado
                  w = [x,y]
        t :  tempo
        p :  vetor de parâmetros conhecidos
                  p = [m,k,c]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, k3, c1, c2, c3 = p
    f1 = 0.
    f2 = 0.
    # Cria a função f = (x',y'):
    sistema = [y1,
               (f1 - (c1+c2)*y1 + c2*y2 - (k1+k2)*x1 + k2*x2 ) / m1,
               y2,
               (f2 - c2*y1 -(c2+c3)*y2 + k2*x1 - (k2+k3)*x2 ) / m2]
    return sistema<jupyter_output><empty_output><jupyter_text>Agora podemos escrever a parte do programa para efetuar a integração do sistema de EDO:<jupyter_code>import numpy as np       # Importa o módulo numpy (métodos numéricos)
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import matplotlib.animation as animation


# variáveis dos enunciado

## Parâmetros da suspensão
m1 = 1.             # Massa
m2 = 1.

k1 = 00.             # Coeficientes de mola
k2 = 50.
k3 = 0.

c1 = 0.              # Coeficientes de amortecimento
c2 = 1.
c3 = 0.


x10 = 1.0             # Condição inicial para x1
y10 = 0.0             # Condição inicial para y1 (velocidade inicial de m1)
x20 = 0.0             # Condição inicial para x2
y20 = 0.0             # Condição inicial para y2 (velocidade inicial de m2)

# Parâmetros e critérios de parada do algoritmo de integração ODE
abserr = 1.0e-8        # Erro absoluto
relerr = 1.0e-6        # Erro relativo
stoptime = 10.0        # Tempo final de integração
numpoints = 501       # Número de pontos de integração

# Cria o vetor com os valores dos tempos de integração para utiliza no plot
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
dt = t[1]-t[0]

# Agrupando os parâmetros e condições iniciais par ao algoritmo:
p = [m1, m2, k1, k2, k3, c1, c2, c3]
w0 = [x10, y10, x20, y20]

# Chamada ao algoritmo de integração 
wsol = odeint(sistema, w0, t, args=(p,),atol=abserr, rtol=relerr)

# Separa a solução nos vetores x e y
x1 = wsol[:,0]
y1 = wsol[:,1]
x2 = wsol[:,2]
y2 = wsol[:,3]


fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0,12), ylim=(-0.5,0.5))
ax.set_aspect('equal')
ax.grid()

line, = ax.plot([], [], 'o-', lw=2, markersize=15)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

print('Iniciando animação...')


def animate(i):
   
    thisx = [0, x1[i]+3 , x2[i]+6, 12]
    thisy = [0, 0, 0, 0]
    line.set_data(thisx, thisy)
           
    time_text.set_text(time_template % (i*dt))
    return  line, time_text


ani = animation.FuncAnimation(
    fig, animate, len(wsol), interval=dt*numpoints, blit=True)
plt.show()

ani.save('final2.mp4')
print('Mostrando animação no notebook...')

from IPython.display import HTML
HTML(ani.to_html5_video())<jupyter_output>Iniciando animação...<jupyter_text>![image.png](attachment:image.png)<jupyter_code># Plotando o resultado
plt.figure()         # Define uma figura para plotar
                     # Figura é formada por dois subplots

plt.plot(t,x1, label="x1(t)")
plt.plot(t,x2, label="x2(t)")
plt.legend()
plt.xlabel('t [s]')
plt.ylabel('x1(t), y2(t)')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13274144651364228, 'fraction_numerical': 0.060848852317020356, 'mean_word_length': 3.161261261261261, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '238130', 'license': 'MIT', 'n_tokens_mistral': 2082, 'n_tokens_neox': 1974, 'n_words': 517, 'path': '.ipynb_checkpoints/06 2 DOF_Mov-Animation-checkpoint.ipynb', 'repo_name': 'nuglock/Vibration', 'size': '4628'}