# Relaxation of Mass-Spring System with Pinned Masses

## 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 *)
(1. - l/EuclideanDistance[mass[[i]], mass[[j]]])*(mass[[i]] -
mass[[j]]), {j, spring[[i]]}];
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];
newpotential = potential[mass];

While[Abs[newpotential - oldpotential] > tol,
mass = mass + step*newdir;
olddir = newdir;
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
]
];

• +1 just for a fantastically written question - refreshing!
– ciao
Mar 25, 2015 at 6:37
• @rasher: Thanks! That's great I'm hearing that from you. Mar 25, 2015 at 6:43
• I wrote an answer that should solve the problem by "hand" with LinearSolve, you could also use directly LinearPrograming but possibly the new finite element abilities of NDSolve, in Mathematica 10, may have the means to do this all automatically, but I am still stuck with Mathematica 9 and have no means to play with it :,( . Mar 25, 2015 at 12:28

(*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};


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.

• Thanks for your elegant answer. I checked this answer against my own, it seems on the more complex networks, mine finds lower energies (though much slower!). Apr 11, 2015 at 22:13
• The built-in Minimize function is a "black box" which, it seems to me, is not really able to find the lowest result in every possible kind of problem. Here is another easy way to minimize the energy in a system like yours: make a small random change in the position of one particle, do a test evaluation, and retain the change if it brings about a lower value for the energy, otherwise undo the change. Repeat this until the energy fails to decrease over a suitably large number of iterations. But if you prefer not to use any randomness in finding an answer, there are still several other methods. Apr 12, 2015 at 5:48

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}];