2D spring-mass systems in equilibrium

Vector notation preliminaries

First, we summarize 2D vector notation used in the derivations for the spring
system.
We use a to denote the length of a vector a, a =
_
a
2
x
+ a
2
y
. The outer
product ab
T
of two vectors a and b is a matrix
_
a
x
b
x
a
x
b
y
a
y
b
x
a
y
b
y
_
Dierentiating a scalar function f(p) of a vector argument p = (x, y), with
respect to p means calculating the vector
f
p
=
_
f
x
,
f
y
_
.
If we dierentiate a vector function with respect to a vector we get a matrix
(one can think about this as dierentiating the vector function component by
component, and stacking two resulting row vectors together); if the components
of f (p) are f
x
(x, y) and f
y
(x, y), then
f
p
=
_
fx
x
fx
y
fy
x
fy
y
_
This notation makes linear approximations to vector functions look exactly
the same way as in the scalar case: instead of
f(x) f(x
0
) + f

(x)(x x
0
)
we get
f (p) f (p
0
) +
f
p
(p p
0
)
where the second term is a matrix-vector product. Several useful identities:

p
p
= I
that is, if f
x
(x, y) = x and f
y
(x, y) = y, the matrix f /p is just the
identity matrix.
For a scalar function g, and vector function f , we get the product rule:
gf
p
= g
f
p
+f
g
p
T
where the second term is an outer product. This can be easily veried
applying the denitions of the derivative of the vector function, the usual
scalar derivative product rule, and the outer product denition.
1
Chain rule for two vector functions g(p) and f (q):
f (g(p))
p
=
f
q
g
p
where the right-hand side is the matrix product.
Spring equations
p
i
p
j
p
j
r
ij
Rest configuration
One point moved
F
ij
F
ji
A spring-mass system is a collection of point masses m
i
with positions p
i
con-
nected by springs. Some of the points are xed, some are allowed to move.
Our goal is to nd positions of the moving points for which the total force from
springs acting on each point is zero, in other words, the system is in its rest
(equilibrium) state to which it relaxes if no external forces are applied.
A spring (i, j) connects two points p
i
and p
j
. Each spring has a rest length,
r
ij
. If the distance p
i
p
j
between spring endpoints is r
ij
, there there is
no contribution from the spring to the total force acting on either point. If
the distance is dierent, there is a spring force acting on each endpoint: the
magnitude of the force acting on p
i
, is kp
i
p
j
, and it acts along the spring,
that is, the line connecting spring endpoints; the coecient k is called stiness.
For simplicity, we assume that all springs have the same stiness. In vector
form, the force can be written as
F
ij
= k(p
i
p
j
r
ij
)
p
j
p
i
p
i
p
j

In this expression, p
j
p
i
/p
i
p
j
is just a unit vector pointing along the
spring away from p
i
.
We denote neighbor(i) the set of indices of points connected to point i by a
spring. Then the total force acting on a point i is
F
i
=

jneighbor(i)
F
ij
=

jneighbor(i)
k(p
i
p
j
r
ij
)
p
j
p
i
p
i
p
j

(1)
Let I
f
be the set of indices of points that are xed, and I
m
be the set
of indices of moving points; we assume that all moving points precede in the
list of points all xed points, that is, the moving points have indices I
m
=
{0 . . . N
m
1}, and xed points I
f
= {N
m
. . . N
m
+ N
f
1}. Later we will
discuss how to eliminate this assumption.
We build a vector F of length 2N
m
, which is obtained by concatenating 2D
vectors of forces F
i
for all moving points. Similarly, we concatenate all moving
point positions p
i
to get a vector p.
2
Then the equation we need to solve has the form
F(p) = 0 (2)
One can see that for the equilibrium problem, the actual mass of each point
does not matter, because the forces only depend on spring stinesses (for the
dynamic problem, F
i
= m
i
a
i
, where a
i
is the acceleration of i-th point, the
mass is important).
These equations are nonlinear, so we cannot solve this system using, for
example, LU decomposition.
Solving the system using Newtons method
The Newton method iteration for (2) (See Chapter 9 of the textbook for details)
is
p
k+1
= p
k
J[F](p
k
)
1
F(p
k
)
where J[F](p
k
) is the Jacobian matrix of F evaluated for point positions p
k
.
The Jacobian matrix is the matix of partial derivatives with entry (s, t) equal
to F
s
/p
t
, where s, t = 0 . . . 2N
m
1, are indices in vectors F and p. Note
that this is a more general case of the 2 2 derivative matrix f /p.
It is convenient to view each Newton step as solving a linear system
J[F]p
k+1
= F(p
k
)
where p
k+1
= p
k+1
p
k
.
We already know how to obtain the vector F(p
k
); so the only other element
that we need is the Jacobian matrix.
Assembling the Jacobian matrix. It is convenient to assemble this matrix
from 2 blocks, obtained by dierentiating individual spring forces F
ij
. Instead
of treating p as a at vector [p
x
1
, p
y
1
, p
x
2
, p
y
2
, . . .], we regard it as a vector of points
[p
1
, p
1
, p
2
. . .], and take the same approach to F. Now we build J[F] as a
N
m
N
m
matrix of 2 2 blocks F
i
/p
j
, i, j = 0 . . . N
m
1.
An individual block F
i
/p
j
can be interpreted as the rate of change of the
force acting on point p
i
due to the change in position of point p
j
.
First, we consider the case when i = j. In this case, in the sum F
i
=

l
F
il
,
at most one term depends on p
j
: if there is a spring connecting p
i
with p
j
,
then F
ij
depends on j, and the rest do not, and if there is no spring, no terms
depend on j. So we get, for i = j,
F
i
p
j
=
_
Fij
pj
, if there is a spring (i, j)
0 otherwise.
If i = j, all terms in the sum depend on p
j
= p
i
, so the derivative will be
the sum:
F
i
p
i
=

lneighbor(i)
F
il
p
i
(3)
On the other hand, we observe that F
il
depends only on the dierence d
il
=
p
i
p
l
, so we can write it as a composition of functions F(d
il
(p
i
p
l
)). As
d
il
/p
i
= I and d
il
/p
l
= I, by the chain rule, we can see that F
il
/p
i
=
F
il
/p
l
.
Based on this, we see that we can rewrite the expression (3) as
3
F
i
p
i
=

lneighbor(i)
F
il
p
l
There is a subtlety here related to xed points: some of the points p
l
attached
to the moving point p
i
may be xed; we still need to include terms F
il
/p
l
in
the summation above, although p
l
for a xed point is not changing.
So all blocks in the matrix can be obtained as soon as we know F
ij
/p
j
.
Computing force derivative for an individual spring. To compute this
derivative, we rewrite F
ij
as sum of two terms
F
ij
= k(p
j
p
i
) kr
ij
p
j
p
i
p
i
p
j

The derivative of the rst term with respect to p

j
is kI. The derivative of
the second term requires a bit more eort. We compute it using the chain rule:
denote d = p
j
p
i
. Then we can write the second term as
G
ij
= kr
ij
d
d
= kr
ij
d

d d
The complete derivative is obtained as (G
ij
/d)(d/p
j
). But d/p
j
=
p
j
/p
j
p
j
/p
i
= I 0 = I, so our derivative is simply G
ij
/d. As kr
ij
is
constant, it remains to compute (d/d)/d, which we do using the product
rule:
d
1
d
d
= d
1
d
d
+d
_
d
1
d
_
T
The derivative
(dd)
1/2
d
is (1/2)(d d)
3/2 dd
d
= d
3
d, so we get the
expression
d
1
d
d
=
d
2
I dd
T
d
3
And the nal expression for the derivative of the force F
ij
with respect to
the point position p
j
is
J
ij
=
F
ij
p
j
= kI kr
ij
d
ij

2
I d
ij
d
T
ij
d
ij

3
where d
ij
= p
j
p
i
.
An additional observation we can make about these expression is that J
ij
=
J
ji
, because d
ij
= d
ji
, and sign changes of d cancel in the expressions. This
allows us to compute only J
ij
for j < i, and use it for two subblocks of the
matrix.
By-spring assembly. The total forces F
i
acting on each point can be com-
puted if we iterate over all points and compute the sums (1) for each. This
would require additional data structures for adjacency information, and a loop
for each point. Alternatively, we observe that we can have a loop over springs
(i, j), and add forces due to each spring to the parts of the total force vector
corresponding to points p
i
and p
j
:
for each spring (i, j), i < j
compute F
ij
if i < N
m
add F
ij
to F
i
if j < N
m
add F
ij
to F
j
.
4
Here we take advantage of the fact that F
ij
= F
ji
to compute the force
only once.
Similarly, instead of assembling the matrix J[F] using a double loop over F
s
for rows and p
t
for columns, we can compute F
ij
/p
j
for each spring (i, j),
and add it to the 2 2 block in the matrix corresponding to the pair of points
(i, j). As we know that the coordinates of i-th moving point occupy slots 2i
and 2i + 1 in the vectors F and p, thanks to the assumption that all moving
points precede xed, the block J
ij
corresponds to the matrix subblock in rows
2i, 2i + 1 and columns 2j, 2j + 1.
This leads to the following algorithm that can be implemented with just sev-
eral lines of code, where we use numpy notation for matrix subblock assignment:
Initialize the matrix J to 2N
m
2N
m
zero matrix
foreach spring (i, j), i < j
compute J
ij
=
Fij
pj
if i < N
m
# update the sum for the diagonal block, if i is not fixed
J[2i : 2i + 2, 2i : 2i + 2]+ = J
ij
if i < N
m
and j < N
m
,
# set off-diagonal blocks if neither point is fixed
J[2i : 2i + 2, 2j : 2j + 2] = J
ij
J[2j : 2j + 2, 2i : 2i + 2] = J
ij
An alternative way of handling xed points. If we x and unx points
interactively, like we do in our code, it is convenient not to assume that indices
of all moving points precede indices of xed points: this would require resorting
the points each time.
The idea is to regard all points, including xed, as variables, but modify
the system so that moving points stay in their positions. That is, now the
vector p has all points in it, no longer in any particular order, and its length is
2(N
m
+ N
f
) = 2N. The force vector F is assembled in the same way as p.
Then we zero out the entries in the force vector corresponding to xed points,
so that there is no force acting on them; as a consequence, these do not move.
Now the matrix J[F] is 2N 2N. We start with assembling it as if all points
were moving: we use the algorithm above, but eliminate the check if the points
p
i
and p
j
are moving or not. Finally, we replace the rows and columns of the
matrix corresponding to xed points with corresponding rows and columns of
identity matrix.
One can verify that for this approach the xed points are not modied by
the Newton iteration, and the submatrix for the moving points has exactly the
same entries as in the original algorithm.
Here is an illustration of this approach in the case of three points, p
0
, p
1
, p
2
,
with p
0
and p
2
xed, and two springs (0, 1) and (1, 2):
1. Initially, we set F = [F
01
, F
01
+F
12
, F
12
]. Then, we set forces on xed
points to zero: F = [0, F
01
+ F
12
, 0].
2. To assemble J[F], we rst assemble the matrix as if all points were moving:
_
_
J
01
J
02
J
01
J
02
J
01
J
01
J
12
J
12
J
02
J
12
J
02
J
12
_
_
5
3. Then we set the rows and columns corresponding to xed points 0 and 2
to the rows and columns of the identity matrix of the same size as J[F]:
J[F] =
_
_
I 0 0
0 J
01
J
12
0
0 0 I
_
_
The inverse of this block-diagonal matrix has the same block form, with
diagonal entries inverted:
J[F]
1
=
_
_
I 0 0
0 (J
01
J
12
)
1
0
0 0 I
_
_
When we multiply J[F]
1
by F for the Newton update step, we get J[F]
1
F =
[0, (J
01
J
12
)
1
(F
01
+F
12
), 0], so there is no change in positions of the xed
points, and the correct update, that we would have obtained using the rst
method for matrix assembly for the moving points. This requires maintaining a
larger matrix, but eliminates the need to change point indices when points are
xed and unxed.
6