Many Driven Springs

First pick a number of moving masses:

num = 7

Next decide to use fixed or periodic boundary conditions,
i.e. whether masses 1 and N are attached to fixed walls, or to each other.

FixedBC = True (*Choose True or False here*)

(mfixed[i_,j_] = If[i==j,2.0,If[(i==j+1)||(i==j-1),-1.0,0]];
mperiodic[i_,j_] = If[i==j,2.0,If[(i==j+1)||(i==j-1)||(i==num&&j==1)||(i==1&&j==num),-1.0,0]];
mm = Table[ If[FixedBC,mfixed[i,j],mperiodic[i,j]],{i,num},{j,num}];
MatrixForm[mm])

Ask for the eigenmodes and their frequencies.
Also define an eigenmode browser.

eig=Eigenvectors[mm];
ω0 = Sqrt[Abs[#]]& /@ Eigenvalues[mm];
draweig[j_]:=ListLinePlot[eig[[j]],AxesLabel->ω0[[j]]]
Manipulate[draweig[n],{n,1,num,1}]

Choose a force vector and find the complex coefficient of the j’th mode in terms of generic driving ω.

fchoice=1;  (* choose which type of force *)
f=Which[
fchoice==1,Table[1-2Random[],{num}],(*literally random*)
fchoice==2,f={0,0,0,1,0,0,0},
fchoice==3,f={1,0,0,0,0,0,0},
fchoice==4,f=eig[[2]] (*pick one of the eigenvectors*)
];
f=f/Norm[f]

Decompose in to eigenmode components and divide by ω^2 etc. to get the resulting amplitude
and phase of oscillations, call it q.  Then the power goes like m ω^2 γ |q|^2.

feig = eig.f
q=Table[feig[[j]]/(-ω^2+ω0[[j]]^2 + 2 I ω γ),{j,num}]
qq=-Im[(eig.f).q/.{γ->.02}];
qq = ω^2 γ Norm[q]^2
Plot[qq/.γ->.04,{ω,0,3},PlotRange->{0,3}]

Here are some drawing functions for 3D springs:

spring[v1_,v2_]:=Module[{q,e1,e2,e3,n=10,z={0,0,1}},(q=v2-v1;e1=q/Norm[q];q=z-(e1.z)e1;e2 = q/Norm[q];e3=Cross[e1,e2];{ParametricPlot3D[v1+s(v2-v1)+.1Tanh[s(s-1)2 π n]*( Cos[2π n s]e2+Sin[2π n s]e3),{s,0,1},PlotStyle->{Red,Thick}][[1]]})]
state3D[t_] := Graphics3D[{
{Black,forceGraphics[t]},
{Green,Thick,Table[Line[{{i,0,-1},{i,0,1}}],{i,num}]},
{GrayLevel[.25],Specularity[LightBlue,10],Table[Sphere[pt[i,t],0.2],{i,num}]},
Table[spring[pt[i,t],pt[i+1,t]],{i,0,num}]
},
Boxed->False,PlotRange->{{0,num+1},{-2,2},{-2,2}},ViewPoint->{0,2,1},ImageSize->800]

And making a final numerical choice of parameters, we Animate:

thisω0=.9; thisγ=.02;
forceGraphics[t_]:= Table[Sphere[{j,0,f[[j]]Cos[thisω0 t]},.05],{j,num}]
coords[t_]:=Transpose[eig].Re[q Exp[I ω t]]/.{ω->thisω0,γ->thisγ};
Do[pt[i,t_]={i,0,coords[t][[i]]},{i,1,num}];
If[FixedBC,(pt[0,t_]={0,0,0};pt[num+1,t_]={num+1,0,0};),(pt[0,t_]=Join[{0},Rest[pt[num,t]]];pt[num+1,t_]=Join[{num+1},Rest[pt[1,t]]];)]
Animate[state3D[t],{t,0,30}]