Harmonic Oscillators

Functions

Below we have the basic solver - it takes in a (symmetric) Q matrix, the constant offset and initial positions and velocities, then diagonalizes and solves the problem using (8.38) and (8.39):

Solveit[Qin_, Ain_, X0in_, DX0in_] := Module[{i, Psi0, DPsi0, psiout,B, V, S},
psiout = {};
 {S,V} =N[ Eigensystem[Qin]];
V =Orthogonalize[V];
V = Transpose[V];
Psi0 = Transpose[V].X0in;
 DPsi0 = Transpose[V].DX0in;
B = Transpose[V].Ain;
For[i = 1, i ≤ Length[S], i = i + 1,
If[S[[i]] ≠ 0,
psiout = Join[psiout,{(1/(2 S[[i]]))(E^(-t (Sqrt[S[[i]]]))) ((((-1+(E^(t (Sqrt[S[[i]]])))))^2) B[[i]]+(-1+(E^(2 t (Sqrt[S[[i]]])))) DPsi0[[i]] (Sqrt[S[[i]]])+(1+(E^(2 t (Sqrt[S[[i]]])))) Psi0[[i]] S[[i]])}];
,
psiout = Join[psiout,{(1/2) (t^2) B[[i]]+t DPsi0[[i]]+Psi0[[i]]}];
];
];
V.psiout
]

Below are some helper functions - they place points at locations given in the list X.

DrawConfig[X_] :=Graphics[ Table[{PointSize[.025],GrayLevel[1-jj/Length[X]],Point[{X[[jj]],0}]},{jj,1,Length[X]}]]

DrawConfig3D[X_] :=Graphics3D[ Table[{PointSize[.025],Point[{X[[3 (jj-1)+1]],X[[3 (jj-1)+2]],X[[3 (jj-1)+3]]}]},{jj,1,Length[X]/3}]]

Three Balls and Two Springs

Here we have three balls conected by two springs with the same equilibrium spacing a and spring constants k

k = 1;
m = 1;
a = 1;

We make the matrix and offset vector

Q = {{-k/m, k/m, 0},{k/m, -(k+k)/m, k/m},{0,k/m, -k/m}};
A = {-k a/m, (k - k) a/m, k a/m};

and provide the initial conditions, chosen arbitarily here

X0 = {-1,0,0};
V0 = {0,0,0};

XX = Solveit[Q,A,X0,V0];

Let's check the solution:

Print[Chop[Simplify[D[D[XX,t],t] -Q.XX - A]]," ",Chop[(XX//.t-> 0 ) - X0]," ",Chop[(D[XX,t]//.t-> 0 ) - V0]]

Visualization

Using the DrawConfig function, we can plot the three points and their motion in time.

ListAnimate[Table[Show[DrawConfig[XX //. t-> T],PlotRange-> {{-5,5},{-1,1}}],{T,0,20,.2}]]

Linearized Springs in 2 and 3 Dimensions

The following takes an equilibrium position and a set of connectivities of the form: {{i1,j1},{i2,j2},{i3,j3}, . . . {im, jm}} specifying connections between particle number i1 and j1 for example, and returns the Hessian for a potential with this connectivity and all springs with the same constant k - the matrix is defined in (8.44) and (8.46):

MakeQ[Xbarin_,k_, Connections_] := Module[{Qout,index,i,j,ttsub,rij,a}, 
Qout = Table[0,{ii,1,3 Length[Xbarin]},{jj,1,3 Length[Xbarin]}];
For[index = 1, index ≤ Length[Connections], index = index+1,
i = Connections[[index,1]];
j = Connections[[index,2]];
rij = Xbarin[[i]] - Xbarin[[j]];
a = Sqrt[rij.rij];
ttsub = Table[rij[[jj]] rij[[kk]],{jj,1,3},{kk,1,3}]/a^2;
Qout[[ 3 (i-1)+1;; 3 (i-1)+1+2, 3 (i-1)+1;; 3 (i-1)+1+2]] += ttsub;
Qout[[ 3 (j-1)+1;; 3 (j-1)+1+2, 3 (j-1)+1;; 3 (j-1)+1+2]] += ttsub;
Qout[[ 3 (i-1)+1;; 3 (i-1)+1+2, 3 (j-1)+1;; 3 (j-1)+1+2]] += -ttsub;
Qout[[ 3 (j-1)+1;; 3 (j-1)+1+2, 3 (i-1)+1;; 3 (i-1)+1+2]] += -ttsub;
];
Qout =  -k Qout;
Return[Qout];
]

Small Example

Start with four balls, and nearest neighbor connections as shown below

Ntot =4;
k = 1; 
a = 1;
Xbar = Partition[Flatten[Table[{ii, jj,0},{ii,1,2},{jj,1,2}]],3];
Xbar0 = {{0,0,0},{1,0,0},{0,1,0},{1,1,0}};
Connections = {{1,2},{1,3},{2,4},{3,4}};
Show[Graphics3D[Table[Text[ToString[jj],Xbar[[jj]]],{jj,1,Length[Xbar]}]]]

Make the matrix, and right-hand side vector from (8.46):

Q = MakeQ[Xbar,k, Connections];
A = -Q.Flatten[Xbar];
X0 = Flatten[Xbar] + Table[RandomReal[{-.1,.1}],{ii,1,3 Ntot}];
V0 = Table[{0,0,0},{ii,1,Ntot}];

Solve, check the solution, and visualize it

XX= Solveit[Q,A,X0, Flatten[V0]];

Print[Chop[Simplify[D[D[XX,t],t] -Q.XX - A]]," ",Chop[(XX//.t-> 0 ) - X0]," ",Chop[(D[XX,t]//.t-> 0 ) - Flatten[V0]]]

ListAnimate[Table[Show[DrawConfig3D[XX //. t-> T],PlotRange-> {{0,3},{0,3},{-1,1}}],{T,0,20,.1}]]

Larger Example

We can do the same thing for a larger set of balls.  The connectivity here describes a linear chain that goes in order from ball 1 to ball 25.  The equilibrium configuration is a square sheet:

Ntot =25;
k = 1; 
a = 1;
Xbar = Partition[Flatten[Table[{ii, jj,0},{ii,1,5},{jj,1,5}]],3];
Connections = Table[{i,i+1},{i,1,Ntot-1}];
Show[Graphics3D[Table[Text[ToString[jj],Xbar[[jj]]],{jj,1,Length[Xbar]}]]]

We’ll take the initial configuration to be the equilibrium configuration (the sheet, shown above) with some small random perturbation.  The configuration starts from rest.

Q = N[MakeQ[Xbar,k, Connections]];
A = -N[Q.Flatten[Xbar]];
X0 = N[Flatten[Xbar]] + Flatten[Table[{RandomReal[{-.1,.1}],RandomReal[{-.1,.1}],RandomReal[{-.1,.1}]},{ii,1,Ntot}]];
V0 = Table[{0,0,0},{ii,1,Ntot}];

XX = Solveit[Q,A,X0, Flatten[V0]];

ListAnimate[Table[Show[DrawConfig3D[XX //. t-> T], PlotRange-> {{-1,6},{-1,6},{-1,1}}],{T,0,40,.4}]]

Spherical Example

For the surface of the sphere, we’ll use connectivitiy defined by distance in equilibrium - given some cutoff and a set of particle locations, the function below gives a connectivity list that consists of all particles within cutoff of each other

MakeCons[inX_, cutoff_] := Module[{rlen,rij, conlist, index,jndex},
conlist = {};
For[index = 1, index ≤ Length[inX], index=index+1,
For[jndex = index+1, jndex ≤ Length[inX], jndex =jndex+1,
rij = inX[[index]] - inX[[jndex]];
rlen = Sqrt[rij.rij];
If[rlen ≤ cutoff,
conlist = Join[conlist,{ { index, jndex}}];
];
];
];
Return[conlist];
]

A specialized ViewParticles function for the spherical distribution of interest:

ViewParticles[inX_,cons_] := Module[{out,pts,lines},
pts = Table[Point[inX[[j]]],{j,1,Length[inX]}];
lines = Table[Line[{inX[[cons[[j,1]] ]],inX[[cons[[j,2]] ]]}],{j,1,Length[cons]}];
out = Show[Graphics3D[{Join[{PointSize[0.01]},pts]}],PlotRange->{{-1.5,1.5},{-1.5,1.5},{-1.5,1.5}}];
Return[Show[out,Graphics3D[lines]]];
]

Set up a spherical shell of points, including the top and bottom (carefully):

r = 1.0;
NT = 10;
NP = 20;
dT = N[Pi/NT];
dP = N[2Pi/NP];
Xbar = Flatten[Table[{r Sin[j dT] Cos[k dP], r Sin[j dT] Sin[k dP], r Cos[j dT]},{j,1,NT-1},{k,1,NP}]];
Xbar = Partition[Xbar,3];
Xbar = Join[Xbar,{{0,0,r}}];
Xbar = Join[Xbar,{{0,0,-r}}];cons = MakeCons[Xbar, 0.5];

We’ll take Xbar to be the equilibrium configuration as well as the initial configuration - for initial velocities, we’ll give a slight kick to the top mass:

X0 = Xbar;
V0 = Table[{0.,0.,0.},{j,1,Length[Xbar]}];
V0[[Length[V0] - 1]] = {0,0,-.1};

ViewParticles[Xbar,cons]

Make the matrix and vector:

Q = N[MakeQ[Xbar,1.0,cons]];
A = -N[Q.Flatten[Xbar]];

Solve and view

XX = Solveit[Q,A,Flatten[X0], Flatten[V0]];

XX = Partition[XX,3];
mov = Table[ViewParticles[Chop[XX/. t-> T,10^(-9)],cons],{T,0.,10.,.5}];

ListAnimate[mov]