Figure 8.8:

A hanging chain at rest. See Exercise~\ref {ex:hanging_chain}\ref {sub:hc}.

Code for Figure 8.8

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
# [makes] hanging_chain_skeleton.mat
# Solve the hanging chain QP (unconstrained problem, without ground constraints).
# N mass points at positions (y_i, z_i) minimise total potential energy subject
# to fixed endpoints.
import numpy as np
from scipy.optimize import minimize
from scipy.io import savemat

N   = 40
m_i = 40.0 / N
D_i = 70.0 * N
g0  = 9.81

# Fixed endpoints.
y_first, z_first = -2.0, 1.0
y_last,  z_last  =  2.0, 1.0

# Optimise over the 2*(N-2) free interior coordinates interleaved as
# [y_1, z_1, y_2, z_2, ..., y_{N-2}, z_{N-2}].
def get_yz(x_free):
    y = np.empty(N)
    z = np.empty(N)
    y[0],  z[0]  = y_first, z_first
    y[-1], z[-1] = y_last,  z_last
    y[1:-1] = x_free[0::2]
    z[1:-1] = x_free[1::2]
    return y, z

def objective(x_free):
    y, z = get_yz(x_free)
    spring = D_i/2 * np.sum((np.diff(y)**2 + np.diff(z)**2))
    gravity = g0 * m_i * np.sum(z)
    return spring + gravity

x0 = np.zeros(2*(N - 2))
result = minimize(objective, x0, method='L-BFGS-B',
                  options={'ftol': 1e-12, 'gtol': 1e-10})

Y0, Z0 = get_yz(result.x)

savemat('hanging_chain_skeleton.mat', {'Y0': Y0, 'Z0': Z0})