Npts=30;(*Number of oscillators*)
theta=Pi/(Npts+1);
xpos=Table[j/(Npts+1),{j,1,Npts}];
w=Table[2Sin[theta* j/2],{j,1,Npts}];

(*Normal modes*)

nmode=3;Show[{Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],0.25*Sin[theta nmode j]}],{j,1,Npts}]},PlotRange->{{-0.1,1.1},{-0.5,0.5}}],Plot[{0.25*Sin[Pi nmode x],0},{x,0,1},PlotStyle-> {{Thick,Dashed},{Thick,Dotted}}],Graphics[{Thick,Line[{{0,-0.05},{0,0.05}}],Line[{{1,-0.05},{1,0.05}}]}]}]

y0=0.1;
A=2y0/(Npts+1);
mcenter=0.5*Npts;(*Gaussian center*)
s=2*(0.05*Npts)^2;
aa=A Table[Sum[Exp[-(j-mcenter)^2/s]Sin[theta*n*j],{j,1,Npts}],{n,1,Npts}];

Animate[
Show[{Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[aa[[nmode]]*Cos[w[[nmode]]t]*Sin[theta*nmode*j],{nmode,1,Npts}]}],{j,1,Npts}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}],Graphics[{Thick,Line[{{0,-0.05},{0,0.05}}],Line[{{1,-0.05},{1,0.05}}]}],
Graphics[{Dashed,Line[{{0,0},{1,0}}]}]
}],
{t,0,4Npts},AnimationRunning->False]

(*Caso amortecido*)
f=0.5;(*f<2 so all modes are under-damped*)
gamma=f*Sin[theta/2];

ww=Table[Sqrt[(2Sin[theta* j/2])^2-gamma^2],{j,1,Npts}];

y0=0.1;
A=2y0/(Npts+1);
mcenter=0.5*Npts;(*Gaussian center*)
s=2*(0.05*Npts)^2;
aa=A Table[Sum[Exp[-(j-mcenter)^2/s]Sin[theta*n*j],{j,1,Npts}],{n,1,Npts}];

Animate[
Show[{Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[aa[[nmode]]*(Cos[ww[[nmode]]t]+gamma Sin[ww[[nmode]]t]/ww[[nmode]])*Exp[-gamma t]*Sin[theta*nmode*j],{nmode,1,Npts}]}],{j,1,Npts}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}],Graphics[{Thick,Line[{{0,-0.05},{0,0.05}}],Line[{{1,-0.05},{1,0.05}}]}],
Graphics[{Dashed,Line[{{0,0},{1,0}}]}]
}],
{t,0,8Npts},AnimationRunning->False]

(*Caso forçado*)

LL=7;(*Position of the forced oscillator*)
Omega=2Sin[theta*2/2];
(*Omega=(2Sin[theta*6/2]+2Sin[theta*7/2])/2;*)
(*Omega=1.2w[[Npts]];*)
FF=0.001*Table[Sin[theta*LL*j],{j,1,Npts}];
twoGW=2gamma*Omega;
AA=-twoGW*Table[FF[[j]]/((Omega^2-w[[j]]^2)^2+(twoGW)^2),{j,1,Npts}];
BB=Table[(w[[j]]^2-Omega^2)*FF[[j]]/((Omega^2-w[[j]]^2)^2+(twoGW)^2),{j,1,Npts}];
bbb=Table[(-gamma*AA[[j]]-Omega*BB[[j]])/ww[[j]],{j,1,Npts}];

Animate[Show[{
If[LL==1,
Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*j],{nmode,2,Npts}]}],{j,LL+1,Npts}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}],
{Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*j],{nmode,1,Npts}]}],{j,1,LL-1}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}],Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*j],{nmode,1,Npts}]}],{j,LL+1,Npts}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}]
}],Graphics[{PointSize[Large],Red,Point[{xpos[[LL]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*LL],{nmode,1,Npts}]}]}],Graphics[{Thick,Line[{{0,-0.05},{0,0.05}}],Line[{{1,-0.05},{1,0.05}}]}],
Graphics[{Dashed,Line[{{0,0},{1,0}}]}]}],
{t,0,64Npts},AnimationRunning->False,DisplayAllSteps->True]

(*Caso forçado 2 - 2 osciladores forçados*)

L1=7;(*Position of the forced oscillator*)
L2=24;(*Position of the forced oscillator*)
eta=0;
Omega=2Sin[theta*2/2];
(*Omega=(2Sin[theta*6/2]+2Sin[theta*7/2])/2;*)
(*Omega=1.2w[[Npts]];*)
f0=0.001;
FFs=f0*Table[Sin[theta*L1*j]+Cos[eta]*Sin[theta*L2*j],{j,1,Npts}];
FFc=f0*Table[Sin[eta]*Sin[theta*L2*j],{j,1,Npts}];
twoGW=2gamma*Omega;
AA=Table[((w[[j]]^2-Omega^2)*FFc[[j]]-twoGW*FFs[[j]])/((Omega^2-w[[j]]^2)^2+(twoGW)^2),{j,1,Npts}];
BB=Table[((w[[j]]^2-Omega^2)*FFs[[j]]+twoGW*FFc[[j]])/((Omega^2-w[[j]]^2)^2+(twoGW)^2),{j,1,Npts}];
bbb=Table[(-gamma*AA[[j]]-Omega*BB[[j]])/ww[[j]],{j,1,Npts}];

Animate[Show[{Graphics[{PointSize[Large],Black,Table[Point[{xpos[[j]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*j],{nmode,1,Npts}]}],{j,1,Npts}]},PlotRange->{{-0.1,1.1},{-1.5y0,1.5y0}}],Graphics[{PointSize[Large],Red,Point[{{xpos[[L1]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*L1],{nmode,1,Npts}]},{xpos[[L2]],Sum[((-AA[[nmode]]*Cos[ww[[nmode]]t]+bbb[[nmode]] Sin[ww[[nmode]]t])*Exp[-gamma t]+AA[[nmode]]Cos[Omega t]+BB[[nmode]]Sin[Omega t])*Sin[theta*nmode*L2],{nmode,1,Npts}]}}]}],Graphics[{Thick,Line[{{0,-0.05},{0,0.05}}],Line[{{1,-0.05},{1,0.05}}]}],
Graphics[{Dashed,Line[{{0,0},{1,0}}]}]}],
{t,0,64Npts},AnimationRunning->False,DisplayAllSteps->True]

(*2D*)

Lx=1;
Ly=1;
f[x_,y_,n_,m_]:=Sin[Pi n x/Lx]Sin[Pi m y/Ly]

DensityPlot[f[x,y,3,2],{x,0,Lx},{y,0,Ly}]

Plot3D[f[x,y,3,4],{x,0,Lx},{y,0,Ly}]

w=Table[Pi Sqrt[(n/Lx)^2+(m/Ly)^2],{n,1,50},{m,1,50}];
b=Table[(( Sin[((n π)/4)] Sin[((n π)/2)])/(n π))*(( Sin[((m π)/4)] Sin[((m π)/2)])/(m π))/w[[n,m]],{n,1,50},{m,1,50}];

psi[x_,y_,t_]:=Sum[b[[n,m]]Sin[w[[n,m]]t]Sin[Pi n x/Lx]Sin[Pi m y/Lx],{n,1,20,2},{m,1,20,2}]

Animate[Plot3D[psi[x,y,t],{x,0,Lx},{y,0,Ly},Axes->False,Boxed->False,PlotRange->{-0.02,0.02},Mesh->None],
{t,0,Pi Sqrt[2]},AnimationRunning->False,DisplayAllSteps->True]

Animate[DensityPlot[psi[x,y,t],{x,0,Lx},{y,0,Ly},Axes->False,PlotRange->{-0.02,0.02},Mesh->None],
{t,0,Pi Sqrt[2]},AnimationRunning->False,DisplayAllSteps->True]