Relaxation of Mass-Spring System with Pinned Masses

Question

Problem

Assume that four masses are connected to each other by springs while two of them are pinned down at their initial positions. Find the position of masses when system is at its minimum (potential) energy.

I've got this idea after reading this article. In a mass-spring system, the energy is equal to:

$$ U = \frac{1}{4} k \sum_{ij} (l_{ij}-l_{0})^2 $$

where each spring is counted twice. $l_{ij}$ is the distance of two masses $i$ and $j$ and $l_{0}$ is the spring rest length. My initial try was to use steepest descent to relax the system which was not the best option, so I switched to conjugate gradient.

The gradient on mass $i$ is:

$$ \vec{g}_{i} = \frac{k}{2} \sum_{j} (1-\frac{l_{0}}{l_{ij}}) (\vec{r}_i-\vec{r}_{j}) $$

where $j$ runs over only on masses connected to $i$. Obviously if one pins down all masses, the problem has no solution unless the system is already in the minimum potential. In any other case, the solution is a square (if we assume that masses are hard disks!). Although the above potential allows masses to overlap.

Code

To be specific, consider the following system:

 SeedRandom[30];
 (* initial coordinates of masses *)
 mass = RandomReal[{0, 3}, {4, 2}]; 
 (* mass 1 is connected to {2,3} and so on*)
 spring = {{1, 2}, {1, 3}, {2, 4}, {3, 4}}; 
 (* spring constant *)
 k = 2.;
 (* spring rest length*)
 l = 1.;
 (* step size *)
 step = 0.02;
 (* tolerance *)
 tol = 10^-10.;
 (* pinned masses *)
 pinned = {2,3};

Then I calculate gradient and energy for the system:

(* gradient for mass i *) 
 grad[i_] := Sum[(0.5*k)*
 (1. - l/EuclideanDistance[mass[[i]], mass[[j]]])*(mass[[i]] - 
  mass[[j]]), {j, spring[[i]]}];
(* gradient vector*)
 gradient[mass_] := Flatten[Last@
Reap@Do[If[MemberQ[pinned, i], Sow[{0., 0.}], Sow[grad[i]]], {i, 
   1, Length@mass}], 1];
(* energy of the system*) 
 potential[mass_] := (0.25*k)*Sum[Total[(EuclideanDistance[mass[[i]], mass[[#]]] - l)^2 & /@ spring[[i]]], {i, 1, Length@mass}];

Then by using conjugate gradient, I try to find the minimum potential:

oldgradient = gradient[mass];
oldpotential = potential[mass];
mass = mass - step*gradient[mass];
olddir = newgradient = gradient[mass];
newpotential = potential[mass];

While[Abs[newpotential - oldpotential] > tol,
  gamma = ((Norm@Flatten@newgradient)/(Norm@Flatten@oldgradient));
  newdir = -newgradient + gamma*olddir;
  mass = mass + step*newdir;
  olddir = newdir;
  oldgradient = newgradient;
  newgradient = gradient[mass];
  oldpotential = newpotential;
  newpotential = potential[mass];
];

And finally, system reaches to this point:

which shows 2 and 3 are held fixed.

Questions

1. Do you know of any internal method which does the same job in Matematica?
2. While typing the question, I found "spring embedding algorithm" in Mathematica documentation, under Graph Drawing Algorithm. This is certainly related, so can we find vertex positions using this method?
3. To speed up the code, any help would be appreciated. I tried to write potential and gradient as compiled functions, but they don't seem faster.

Although I am aware that calculating energy of each spring twice is not efficient.

cpotential = With[{w = spring},
 Compile[{{z, _Real, 2}},
    0.25*k*
    Sum[Total[(EuclideanDistance[z[[i]], z[[#]]] - l)^2 & /@ 
    w[[i]]], {i, Length@z}],
    RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed", 
    Parallelization -> True
    ]
];

Answer 1

Here is your initialization code:

(*SeedRandom[30];*)
(*initial coordinates of masses*)

mass = RandomReal[{0, 3}, {4, 2}];
(*mass 1 is connected to {2,3} and so on*)

spring = {{2, 3}, {1, 4}, {1, 4}, {2, 3}};
(*spring constant*)
k = 2.;
(*spring rest length*)
l = 1.;
(*step size*)
step = 0.02;
(*tolerance*)
tol = 10^-10.;
(*pinned masses*)
pinned = {2, 3};

And your potential function:

potential[mass_] := (0.25*k)*
   Sum[Total[(EuclideanDistance[mass[[i]], mass[[#]]] - l)^2 & /@ 
      spring[[i]]], {i, 1, Length@mass}];

Here is a function to display the current position of your masses:

showPositions[positions_] := 
 Graphics[
  {
   AbsolutePointSize[7], Red, Point /@ positions, Black, 
   Table[Style[Text[ToString[j], positions[[j]]], 20], {j, 4}]
  }
 ]

(* Display your initial positions *)

showPositions[mass]

There are many ways to minimize the energy. Here is one:

best = Minimize[
   potential[{{x1, y1}, mass[[2]], mass[[3]], {x4, y4}}], 
   {x1, y1, x4, y4}
];

massOptimized = 
    {{x1, y1}, mass[[2]], mass[[3]], {x4, y4}} /. best[[2]];

showPositions[massOptimized]

The result looks like a parallelogram in which the edge lengths are all 1.

Answer 2

This method is based on the original article. First we need to convert list of spring connectivity (spring) into a graph:

g = Graph[Apply[UndirectedEdge, spring, 1], VertexCoordinates -> mass];

Then Kirchhoff (Laplacian/admittance/discrete Laplacian) matrix should be found using KirchhoffMatrix:

m = KirchhoffMatrix[g] // Normal;

We also need to incorporate pinned coordinates into m by replacing the corresponding rows to zero except the diagonal element:

f[i_] := Array[If[#1 == i, 1, 0] &, Last@Dimensions[m]];
m = m /. (m[[#]] -> f[#] & /@ pinned)

Finally the problem is reduced to solve a matrix equation ($m\mathbf{r}=\mathbf{b}$) in $x$ and $y$ direction:

bx = If[MemberQ[pinned, #], mass[[#, 1]], 0] & /@Range[First@Dimensions[m]];
by = If[MemberQ[pinned, #], mass[[#, 2]], 0] & /@Range[First@Dimensions[m]];
xrelaxed = LinearSolve[m, bx];
yrelaxed = LinearSolve[m, by];

And relaxed positions are simply:

relaxedposition = Thread[{xrelaxed, yrelaxed}];