In [1]: 
%matplotlib inline
In [2]: 
from pylab import *

Un poco de teoría

Las ecuaciones en formato matricial son: \[\overline{\ddot{\Psi}}+\Gamma\overline{\dot{\Psi}}+\mathbb{M}\overline{\Psi}=\overline{F}/m\] En la base de modos normales el sistema se simplifica a \(N\) ecuaciones independientes para cada modo \(i\) de la forma: \[\ddot{\Psi}_i+\Gamma\dot{\Psi}_i+\omega_{0i}^2\Psi_i=F_i/m\] donde \(\omega_{0i}^2\) son los autovalores de \(\mathbb{M}\) y los desplazamientos y fuerzas en los modos se obtienen de \(F=\mathbb{D}^t\overline{F}\) y \(\Psi=\mathbb{D}^t\overline{\Psi}\). La matriz \(\mathbb{D}\) se forma acomodando los autovectores de \(\mathbb{M}\) en columnas.

Solución estacionaria

Para una fuerza proporcional al \(\cos(\Omega t)\), la solución estacionaria para cada nodo esta dada por: \[\Psi_i^{est}=A_i \sin(\Omega t) + B_i \cos(\Omega t)\] con las amplitudes dadas por las relaciones: \[A_i = \frac{F_i}{m}\frac{\Omega\Gamma}{(\omega_{0i}^2-\Omega^2)^2+\Omega^2\Gamma^2}\] \[B_i = \frac{F_i}{m}\frac{\omega_{0i}^2 -\Omega^2}{(\omega_{0i}^2-\Omega^2)^2+\Omega^2\Gamma^2}\] La solución estacionaria para cada partícula se obtiene invirtiendo el cambio de base \(\overline{\Psi}^{est} = \mathbb{D}\Psi^{est}\).

Parámetros del problema

In [3]: 
g = 10                 # aceleración de la gravedad
m = 1                  # masa de los péndulos
l_largos = 1           # largo de los primeros péndulos
l_cortos = 0.5         # idem de los del final de la cadena
k = 5                  # constante de los resortes
Gama = 0.1             # constante de amortiguamiento
N_largos = 30          # número de péndulos largos
N_cortos = 200         # número de péndulos cortos
F0_sobre_m = 4         # fuerza sobre la primera masa
Omega = 3              # frecuencia angular de forzado

Importamos algunas herramientas

In [4]: 
from scipy import *
from scipy.linalg import eigh
from scipy.sparse import spdiags, hstack, vstack, bmat, dia_matrix

Esta función construye la matriz M

In [5]: 
def get_M(N1,l1,N2,l2):
    Nm = N1+N2
    diag0 = zeros(Nm); diag1 = zeros(Nm)
    diag0[:N1] = g/l1+2*k/m          # Elementos de la diagonal, largos
    diag0[N1:] = g/l2+2*k/m          # Elementos de la diagonal, cortos
    diag1[1:] = -k/m                 # Elementos extradiagonales

    Mi = spdiags([diag0,diag1],[0,1],Nm,Nm)
    return Mi

Esta otra función calcula los coeficientes de la solución estacionaria

In [6]: 
def get_sol_estacionaria(w0,Fm,Omega):
    Aes = zeros(w0.size)
    Bes = zeros(w0.size)

    Aes =         Omega*Gama*Fm/((w0**2-Omega**2)**2 + Omega**2*Gama**2)
    Bes = (w0**2 - Omega**2)*Fm/((w0**2-Omega**2)**2 + Omega**2*Gama**2)

    return Aes, Bes

Obtengamos las frecuencias naturales y sus autovectores

In [7]: 
M = get_M(N_largos,l_largos,N_cortos,l_cortos).toarray()
w0_cuadrado, Dmatrix = eigh(M,lower=False)

Ahora calculemos la fuerza sobre los modos

In [8]: 
F = zeros(N_largos+N_cortos); F[0] = F0_sobre_m
F_en_modos = dot(transpose(Dmatrix),F)

Finalmente grafiquemos para algunas frecuencias de forzado las amplitudes

In [9]: 
for i in range(13):
    Omega = 3 + i*3 

    A_en_modos, B_en_modos = get_sol_estacionaria(w0_cuadrado,F_en_modos,Omega)

    A = dot(Dmatrix,A_en_modos)

    B = dot(Dmatrix,B_en_modos)

    figure()
    plot(A,'ok-',label='Amplitud absorbente')
    plot(B,'sr-',label='Amplitud elastica')
    plot((N_largos,N_largos),(1.2*min(A.min(),B.min()),1.2*max(A.max(),B.max())),'k-')
    grid()
    axis([0,40,1.2*min(A.min(),B.min()),1.2*max(A.max(),B.max())])
    xlabel('Masa')
    ylabel('Amplitud')
    legend(loc='lower right')
    w0_min = str(round(w0_cuadrado.min(),2)) 
    w0_max = str(round(w0_cuadrado.max(),2)) 
    title(r'$\Omega = $'+str(Omega)+', '+r'$\omega_0^{min} = $'+w0_min+', '+r'$\omega_0^{max} = $'+w0_max)

Puede explicar lo que observa en los gráficos?

Esperemos que si! :)