# From hanging chains to Ead's bridge

Modelling Question

• What's the shape created by a chain hung between two pins?
• How straight can you make the chain by pulling on one end?

## History of the problem

One of the greatest mathematical puzzles of the 1600's was to find the equation for the shape of a chain hanging between two fixed points. Like many puzzles, nobody really recognized at first that it was a puzzle. We can trace the roots back to one of Leonardo da Vinci's notebooks. Leonardo sketched the shape of a hanging chain, and erroneously asserted that if the supports are not of equal height, then the lowest point of the chain "...will touch the earth nearer to the lower support than to the higher, and in the proportion that the height of the lesser goes into the greater." From there, a succession of natural philosophers, including Isaac Beeckman, Mersenne, Galileo, Huygens, Pardies, and Hooke ruminated on the shape of the hanging chain and the related problem of the suspension bridge.

Fix high up on a wall two nails equally distant from the horizontal ... and from them hang a little thin chain ...; this little chain will bend itself into a parabolic figure. -- Galileo, Discorsi, 1638

And now let this problem be proposed: To find the curve assumed by a loose string hung freely from two fixed points. I assume also that the string is a line which is easily flexible in all its parts -- James Bernoulli, Acta Eruditorum, May, 1690

There appears to have been some confusion about the distinction between the bridge and chain problems, as Galileo repeated the common belief of his time that the shape of a hanging curve was a parabola. This is approximately true in some cases, but not in general. The now-classic solution arises from a problem proposed by Jakob Bernoulli from 1691. Robert Hooke figured out that if you flip a hanging chain over, it will be the "perfect" shape for an archway, as the forces from the weight of the stone itself flow along the curve of the arch. Clifford Truesdell wrote a very interesting history of the subject.

## Empirical observations

Refs: 1, 2, 3, 4

### Dimensional analysis

Suppose we need to hang a rope between two telephone poles and we want to make sure the rope won't come down too close to garage roof. The more slack is in the rope, the farther down it will hang, but we cannot tell how-far down because the rope hangs with a curve. But we can use dimensional analysis to find a relation between the lowest point of a rope hanging between two anchor points of equal height. Call the depth $$d$$ the height from the top of the anchor points to the lowest point on the rope. Well, $$d$$ might depend on the length of the rope $$s$$, the width $$w$$ between the anchors, the density of the rope $$\rho$$, and the acceleration of gravity $$g$$ through a general function relationship $f( d, w, s, \rho, g) = 0.$ The dimensions of a variables will be denoted square brackets [], and the fundamental units will be represented by letters (L)ength, (M)ass, and (T)ime. So,

• $$[d] = L^1 M^0 T^0$$,
• $$[w] = L^1 M^0 T^0$$,
• $$[s] = L^1 M^0 T^0$$,
• $$[\rho] = L^{-1} M^1 T^0$$,
• $$[g] = L^{1} M^0 T^{-2}$$,

Now, let's try to analyze this equation. It is a very general formula – it isn't even solved for one of the variables – but the variables do have units, and that helps. Maybe we can convert the units of these variables to something convenient. We can use conversion factors $$\lambda _ M$$ for units of mass, $$\lambda _ L$$ for units of length, and $$\lambda _ T$$ for units of time. Then, whatever conversion factors we pick, $f( d \lambda _ L, w \lambda _ L, s \lambda _ L, \rho \lambda _ M \lambda _ L^{-1}, g \lambda _ L \lambda _ T^{-2}) = 0$ Close inspection shows that time and mass each only appear once in this equation. With some thought, taking $$\lambda _ L = 1/s$$ , $$\lambda _ M = {1}/{\rho s}$$ , $$\lambda _ T = \sqrt{g/s}$$, will lead to $f\left( \frac{d}{s}, \frac{w}{s}, 1, 1, 1 \right) = 0.$ So, the rope's depth $$d$$ will only depend on the rope's length $$s$$ and the distance between the anchors $$w$$. The density $$\rho$$ of the rope and the strength of gravity $$g$$ have no effect. And this makes some sense -- heavy chains and pieces of thread all seem to hang the same way in common experience. The other feature here is that the aspect ratios $$d/s$$ and $$w/s$$ are dimensionless numbers -- for a given physical problem, these same aspect ratios will give the same number, no matter what units we choose for our measurements!

Solving for $$d$$, and using a new function $$\phi()$$ to express the unknown dependence, $d = s \phi\left(\frac{w}{s} \right).$

If the rope's length equal the width between the anchors ($$s=w$$), the rope has to be perfectly straight, so we should have $$d=0$$. On the other hand, when the anchors are right next to each other ($$w=0$$), we expect the depth to be half the length ($$d = s/2$$). These imply $$\phi(0) = 1/2$$, $$\phi(1) = 0$$, and $$\phi()$$ is a monotone decreasing function in between. So, we have determined quite allot about the shape of a hanging rope without doing very much math.

But as a heads-up, you should be warned that we got rather lucky in choosing to normalize our variables by the rope's length rather than the width between anchors. If we had normalized by the width $$w$$, we would arrive at the alternative dimensionless relationship $$d = w \psi\left(\frac{s}{w} \right)$$. While algebraically equivalent to our first dimensionless relationship for the depth $$d$$, this form is indeterminate when $$w=0$$; we would get in that case $$0 \times \psi(\infty)$$. Our normalization by the rope length $$s$$ avoided this issue since for all physically meaningful versions of the problem, the length $$s$$ will always be positive. For new problems, the best choice of normalization factors may not initially be obvious, and it may be worth-while to consider more than one possibility.

## Derivation of an equation for shape

There are several different approaches for trying to come up with an equation for the shape of a hanging chain.

• Balance of infinitesimal forces.
• Golumb provides a derivation based on tension.(slicker, if you understand tension, and applies to other situations, including elastic beams)
• Karman and Biot's derivation based on moments. (slicker, if you understand moments)
• Functional calculus derivation as a constrained minimizer of potential energy. (A heavy-duty tool in Lagrangian mechanics, but with a lot of extra math)

We'll use the first approach, because it works directly, without first having to develop any new mechanical or mathematical concepts beyond Newton's laws. Let us imagine our chain has very fine small links, so it appears smooth and very flexible (perhaps more like a cable than a chain). This is so we don't have to get in arguments about angles associated with individual pairs of links, but we should keep in mind as a possible source of errors.

So, consider a chain of length $$s$$ hanging between two supports a distance $$w$$ apart.
The length, of course, must be creater than the support separation ($$s>w$$) The chain is fully at rest, not moving. At each point along the chain, there is no acceleration, so the forces of gravity and tension in the chain must be perfectly balanced. We can write this balance of forces as

$F _ {left} + F _ {right} + F _ {gravity} = \vec{0}$

So, let's determine these forces. Let $$y(x)$$ be the height of the chain at each location $$x$$, and that the chain has a uniform thickness and density $$\rho$$ everywhere. At each point $$x$$, there is a little piece of chain of length $$2 \mathrm{d} x$$, $$\mathrm{d} x$$ to the left and $$\mathrm{d} x$$ to the right. The forces are tension to the left and right, and gravity pulling down on that part of the chain. The forces from the tension in the chain must point approximately along the chain. Near the point $$x$$, The dominant parts of the tension forces are then $F _ {left} = k _ L(x)\begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix}$ and $F _ {right} = k _ R(x) \begin{pmatrix} 1 \\ \frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix}$ where $$k _ L(x)$$ and $$k _ R(x)$$ are the unknown magnitudes of the forces. Assuming there is a sizable piece of chain to both the left and right of $$x$$, though, we expect $$k _ L(x) \gg 0$$ and $$k _ R(x) \gg 0$$.

Now, let us look at gravity. We know that gravity must be important because without gravity, the chain could hang in arbitrary shapes or flap in the wind. Assume we are looking at a very short pieces of the chain that is nearly straight, then the force from gravity will be $F_{gravity}=\rho g 2 \mathrm{d} s \begin{pmatrix} 0 \\ -1 \end{pmatrix}$ where $$2 \mathrm{d} s$$ is the length of the piece of chain between $$x-\mathrm{d} x$$ and $$x+\mathrm{d} x$$. Using the arc-length derivation from calculus, we can show that $$\mathrm{d} s = \mathrm{d} x \sqrt{1+y'^2}$$. Therefore our force-balance equation gives $\begin{gather} F _ {left} + F _ {right} + F _ {gravity} = \vec{0} \\ k _ L(x)\begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix}+k _ R(x) \begin{pmatrix} 1 \\ \frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix}+\rho g 2 \mathrm{d} s \begin{pmatrix} 0 \\ -1 \end{pmatrix}=\begin{pmatrix}0\\0\end{pmatrix} \end{gather}$

Since the length of rope is very short, the force of gravity is very small relative to the two tension forces. Thus, the main part of the force balance is $k_L(x) \begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix} + k_R(x) \begin{pmatrix} 1 \\ \frac{\mathrm{d} y}{\mathrm{d} x} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$ This implies that we must have $$k_L(x) = k_R(x)$$. Further, this must be the same everywhere between the two ends of the chain, so this constant should not depend on $$x$$. Thus, we can can just call it $$k$$, with the understanding that $$\partial k/\partial x = 0$$.

We know that gravity must be important because without gravity, the chain could hang in arbitrary shapes or flap in the wind. What we have done so far assumes the chain is straight, but in reality it must be curved.

$F _ {left} = k \begin{pmatrix} -1 \\ -y'(x-dx) \end{pmatrix} = k \begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x \end{pmatrix},$

Using the theory of Taylor series, $y(x+\mathrm{d} x) = y(x) + \frac{\mathrm{d} y}{\mathrm{d} x} \mathrm{d} x + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \frac{1}{2} \mathrm{d} x^2 + \ldots \quad \text{and} \quad y(x-\mathrm{d} x) = y(x) - \frac{\mathrm{d} y}{\mathrm{d} x} \mathrm{d} x + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \frac{1}{2} \mathrm{d} x^2 + \ldots$

Substituting, $F _ {left} = k \begin{pmatrix} -1 \\ -y'(x-dx) \end{pmatrix} = k \begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x+\textrm{h.o.t} \end{pmatrix},$ and $F_{right} = k \begin{pmatrix} 1 \\ y'(x+dx) \end{pmatrix} = k \begin{pmatrix} 1 \\ \frac{\mathrm{d} y}{\mathrm{d} x} + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x+\textrm{h.o.t} \end{pmatrix}$ Back to our force balance, $F_{left} + F_{right} + F_{gravity} = k \begin{pmatrix} -1 \\ -\frac{\mathrm{d} y}{\mathrm{d} x} + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x+\textrm{h.o.t} \end{pmatrix} + k \begin{pmatrix} 1 \\ \frac{\mathrm{d} y}{\mathrm{d} x} + \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x+\textrm{h.o.t} \end{pmatrix} + \rho g 2 \mathrm{d} s \begin{pmatrix} 0 \\ -1 \end{pmatrix}$ $= \begin{pmatrix} 0 \\ 2 k \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x - \rho g 2 \mathrm{d} x \sqrt{1+\left(\frac{\mathrm{d} y}{\mathrm{d} x}\right)^2}+\textrm{h.o.t} \end{pmatrix}$ Focusing now on the vertical component, dividing out the $$\mathrm{d} x$$ and letting $$dx\rightarrow 0$$ (so the higher order terms $$\rightarrow 0$$), and rearranging the terms, we find $\begin{gather*} \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} = \frac{\rho g}{k} \sqrt{ 1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2 } \\ \left( \frac{1}{H} \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \right)^2 = 1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2, \end{gather*}$ where $$H = \rho g/k$$. This is the chain equation. Note that we have not determined $$k$$ -- it is an unknown constant still, and its value depends on the length of the chain and how far the chain anchors are separated.

The chain equation may be solved in closed-form, despite being nonlinear. The standard properties of hyperbolic sine and cosine show that $$\cosh^2 \theta = 1 + \sinh^2 \theta$$ and that $$\dfrac{d}{d\theta} \sinh(\theta) = \cosh(\theta)$$. Fro these properties, it is a direct calculation to show $\begin{gather} \frac{1}{H} \cosh( H (x + C) ) + K \end{gather}$ is a solution for any values of $$C$$ and $$K$$. Note that neither constant $$C$$ or $$K$$ effects the shape of the catenary, only the position. They can be determined by using the positions and heights of the chains two supports. We call these conditions boundary conditions. If we have a chain held up by two supports of equal height and separated by distance $$w$$, then we have boundary conditions $$y(-w/2)=y(w/2)=0$$. Solving for these constants we find $$C=0$$, $$K = -\cosh(Hw/2)/H$$, and $y(x) = \frac{ \cosh( H x ) - \cosh(Hw/2) }{H}$ The shape is controlled completely by $$H$$. To find $$H$$, we need to make sure the hanging chain as the same length as it did when we started. Thus, \begin{align*}s &= \int_{-w/2}^{w/2} \sqrt{ 1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2 } \mathrm{d} x \\ &= \int_{-w/2}^{w/2} \frac{1}{H} \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \mathrm{d} x \\ &= \frac{2}{H} \sinh\left(\frac{wH}{2}\right) \end{align*} With a little rearrangement, we have $\frac{s H}{2}= \sinh\left( \frac{w H}{2} \right).$ This is a transcendental equation for $$H$$ in terms of the chain length $$s$$ and does not have an explicit solution, but can be solved graphically and using 1d root finding as long as $$s > w$$. It is helpful to thing of the constant $$H$$ as an implicitly defined function of the two fundamental problem parameters -- the chain length $$s$$ and the support separation $$w$$, $$H(s,w)$$. The function $$H(s,w)$$ will have three useful properties: $$H(w, w) = 0$$, homogeneity $$H(\lambda s, \lambda w) = \frac{1}{\lambda} H(s,w)$$, and $$\lim _ {s\rightarrow \infty} 1/H(s,w) = 0$$.

## Derivation by minimum potential energy

A fancier approach to this problem is based on Aristotle's idea that nature always does what's best. In this case, we can use Lagrangian mechanics to find the shape of the chain that minimizes the potential energy under the constraint that the chain has fixed length.

The overall idea is to write down an equation for the potential energy of the chain and take the derivative of this with respect to changes in chain shape. Since the derivative should vanish at the minimum potential energy, we set the derivative to zero and solve for the shape. That is also exactly how we would have approached the paroblem in multi-variable calculus. The approach will also work here, with a extra elbow grease needed because we are differentiating with respect to a curve rather than a vector.

### Derivation by functional differentiation

(see Wan's book and Lancos's book for reference)

To minimize the potential energy of the hanging chain

$P = \rho g \int _ 0 ^ w y(x) \sqrt{1+\left(\frac{dy}{dx}\right)^2} dx$

subject to the constraint

$s = \int _ 0 ^ w \sqrt{1+\left(\frac{dy}{dx}\right)^2} dx$

We construct the functional with lagrange multiplier L (which is a constant),

$\rho g \int _ 0 ^ w y(x) \sqrt{1+\left(\frac{dy}{dx}\right)^2} dx + L \left(s - \int _ 0 ^ w \sqrt{1+\left(\frac{dy}{dx}\right)^2} dx \right)$ $\int _ 0 ^ w (\rho g y(x)-L) \sqrt{1+\left(\frac{dy}{dx}\right)^2} dx + s L$ From this, we functionally differentiate with respect to $$y(x)$$ to derive the Euler--Lagrange equation \begin{align} 0 = & \frac{\delta}{\delta y} \int _ 0 ^ w (\rho g y(x)-L) \sqrt{1+y'^2} dx + \frac{\delta}{\delta y} (L s ) \\ = & \int _ 0 ^ w \frac{\delta}{\delta y}\left\{(\rho g y(x)-L) \sqrt{1+y'^2} \right\}dx + 0 \\ = & \int _ 0 ^ w \frac{\delta}{\delta y}\left\{(y(x)-\frac{L}{\rho g}) \sqrt{1+y'^2} \right\}dx \\ = & \int _ 0 ^ w \frac{\delta}{\delta y}\left\{(y(x)- L ) \sqrt{1+y'^2} \right\}dx \\ = & \int _ 0 ^ w \delta(x-u) \sqrt{1+y'^2} + (y(x)-L) \frac{1}{2} (1+y'^2)^{-1/2} \left( 2 y' \delta'(x-u) \right) dx \\ = & \sqrt{1+y'^2} - \frac{\partial}{\partial x} \left\{ \frac{(y-L)y'}{\sqrt{1+y'^2}} \right\} \\ = & \sqrt{1+y'^2} - \frac{(y- L) y'' + \left(y'\right)^{4} + \left(y'\right)^{2} }{\left(\left(y'\right)^{2} + 1\right)^{\frac{3}{2}}} \\ 0 = & \frac{(L-y) y'' + y'^2 + 1}{(1+(y')^2)^{3/2}} \end{align} We note that the denominator is always positive (never zero) and bounded above as long as the first derivative is bounded, in which case we can cross-multiply to get $0 = (L-y) y'' + y'^2 + 1.$ This nonlinear second-order equation. It is NOT the same as the equation we derived from first principles -- our original version included the square of the second derivative, while this version does not. But it turns out the equations have the same kind of general solution, $y(x) = L + \frac{1}{H} \cosh( H (x + C) )$ where $$H$$ and $$C$$ are unknown constants constrained by our arc-length condition above.

## Simulations

### Elastic string under gravity

A simple implementation, using a sequence of nodes on the string.

[Show code]
 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 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96  from numpy import array, arange, linspace, zeros from scipy import sqrt, sin, pi from matplotlib.pyplot import figure,show,xlim,ylim,title, savefig from matplotlib import animation from scipy import integrate fig = figure(figsize=(10.,10.), dpi=72) ax = fig.add_subplot(111) def main(): g = 1.0 # strength of gravity f = 0.1 # friction constant # 23 node chain, fix first and last positions N = 23 resting = 0.2/N # resting distance rho = 1./N # node mass k = 4.0*N # spring constant def vfield(X,t): # x,y give position of node i # u,v give velocity of node i x = X[0*N:1*N] y = X[1*N:2*N] u = X[2*N:3*N] v = X[3*N:4*N] # rate of change in position is given by velocities Dx = u Dy = v # For total force on each node, we start with the # external forces Du = -f*u Dv = -f*v - g # boundary conditions Du[0] = 0 Dv[0] = 0 Du[N-1] = 0 Dv[N-1] = 0 for i in range(1,N-1): # internal forces only between nearest neighbors for j in ( i -1, i + 1 ): # do a linear Hooke's law for force # being proportional to the difference between # resting and current length s = x[j] - x[i], y[j] - y[i] distance = sqrt(s[0]**2 + s[1]**2) s = s[0]/distance, s[1]/distance Du[i] += s[0]*k*(distance-resting) Dv[i] += s[1]*k*(distance-resting) # repack all the pieces to be returned result = zeros(4*N) result[0*N:1*N] = Dx result[1*N:2*N] = Dy # normalize forces by node masses to get accelerations result[2*N:3*N] = Du/rho result[3*N:4*N] = Dv/rho return result start = 0 end = 5 dt = 0.01 time=arange(start,end,dt) X = array([0.]*(4*N)) X[0:N] = linspace(0,1,N) X[0:N] = X[0:N] + 0.1*sin(2*pi*X[0:N]) Y = integrate.odeint(vfield,X,time) def animate(i): x = animate.Y[i,0:N] y = animate.Y[i,N:2*N] animate.line.set_data(x,y) animate.title.set_text('Time t=%03.2f'%animate.time[i]) savefig('frame%04d.tif'%i) ## Bind our grid to the identifier X in the animate function's namespace. animate.line, = ax.plot(Y[0,0:N],Y[0,N:(2*N)],'ko-') animate.title = title('Time t=%03.2f'%0) animate.time = time animate.Y = Y xlim(-0.2,1.2) ylim(-1,0.2) # I don't know why it is necessary to create this object. anim = animation.FuncAnimation(fig, animate, len(time), interval=10) show() main() 

Consider using a potential-formulation for force-distance calculations.

Solution of the system of ordinary differential equations describing the motion of an elastic rope

### Truss

Use a graph-model, which provides a more general solution and can accomodate various truss designs and node masses.

[Show code]
 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 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144  from numpy import array, arange, linspace, zeros from scipy import sqrt from matplotlib.pyplot import figure,show,xlim,ylim,title,savefig from matplotlib import animation from scipy import integrate fig = figure(figsize=(10.,10.), dpi=72) ax = fig.add_subplot(111) def gettrussseq(n): # when n = 4, # # 0---1---2---3 # \ / \ / \ / # 4---5---6 # # should give the reordered node sequence # so that we can walk all edges exactly once. # # (4,1,5,2,6,3,2,1,0,4,5,6) # L = [] for i in range(n-1): L.append(i+n) L.append(i+1) L.extend(range(n-2,-1,-1)) L.extend(range(n,2*n-1)) #print L return L def trussreorder(Y): # reorder the ODE solution columns so that # we can draw the entire truss with a single # line and a single animation data update. N = Y.shape[1]/4 n = (N+1)/2 t = gettrussseq(n) M = len(t) Z = zeros((Y.shape[0],M*2)) for (i,j) in zip(range(M),t): Z[:,i] = Y[:,j] Z[:,i+M] = Y[:,j+N] #print n, N, M, Y.shape, Z.shape return Z, M def main(): g = 1.0 # strength of gravity f = 0.2 # friction constant # Create a graph representing a simple truss using a # numbering system like the following. # # 0---1---2---3 # \ / \ / \ / # 4---5---6 # n = 10 # truss length N = 2*n-1 edges = [] edges += [ (i,i+1) for i in range(n-1) ] edges += [ (i,i+1) for i in range(n,2*n-2) ] edges += [ (i,i+n) for i in range(n-1) ] edges += [ (i,i+n-1) for i in range(1,n) ] edges += [ (e[1],e[0]) for e in edges ] edges = list(set(edges)) edict = {} for e in edges: if e[0] not in edict: edict[e[0]] = [] edict[e[0]].append(e[1]) lax = 1/float(n-1) # resting distance rho = 1./n # node mass k = 40.0*n # spring constant def vfield(X,t): x = X[0*N:1*N] y = X[1*N:2*N] xv = X[2*N:3*N] yv = X[3*N:4*N] D = array([0.]*(4*N)) D[0*N:1*N] = xv D[1*N:2*N] = yv for i in range(N): # boundary conditions if i in (0,n-1): continue # external forces D[2*N+i] = -f*xv[i] D[3*N+i] = -f*yv[i] - g # internal forces for j in edict[i]: s = x[j] - x[i], y[j] - y[i] dist = sqrt(s[0]**2 + s[1]**2) s = s[0]/dist, s[1]/dist D[2*N+i] += s[0]*k*(dist-lax) D[3*N+i] += s[1]*k*(dist-lax) # normalize accelerations by node mass D[2*N+i] = D[2*N+i]/rho D[3*N+i] = D[3*N+i]/rho return D start = 0 end = 10 dt = 0.01 time=arange(start,end,dt) # Initial condition X = array([0.]*(4*N)) X[0:n] = linspace(0,1,n) X[n:N] = linspace(0,1,n-1)*(n-2)*lax + lax/2 X[(N+n):(2*N)] = -sqrt(3)/2*lax def animate(i): x = animate.U[i,0:M] y = animate.U[i,M:2*M] animate.line.set_data(x,y) animate.title.set_text('Time t=%03.2f'%animate.time[i]) #savefig('frame%04d.tif'%i) ## Bind our grid to the identifier X in the animate function's namespace. animate.U, M = trussreorder(integrate.odeint(vfield,X,time)) #animate.U = integrate.odeint(vfield,X,time) Y = animate.U animate.line, = ax.plot(Y[0,0:M],Y[0,M:(2*M)],'ko-') animate.title = title('Time t=%03.2f'%0) animate.time = time xlim(-0.2,1.2) ylim(-1,0.2) # why is it necessary to create this object? anim = animation.FuncAnimation(fig, animate, len(time), interval=10) show() main() 

Solution of the system of ordinary differential equations describing the motion of a truss with massless struts and uniform node masses

# Exercises

1. In class, we took a photograph of the arc created by a hanging chain, and asked “What determines this curve’s shape?”. http://www.math.psu.edu/treluga/450/hangingchain.png is photo we had taken. The easiest way to think of a camera image is as a sequences of 3 matrices, where each matrix represents the intensity of one color, red, green, or blue ( sometimes, people will use the word “tensor”). Now, using a little python magic (that you may find useful to read and explore), we can extract data from the image like http://www.math.psu.edu/treluga/450/thresholds.png. http://www.math.psu.edu/treluga/450/extractchaindata.py. , I extracted a data set of points on the curve, which you can read at http://www.math.psu.edu/treluga/450/hw04chaindata.txt

1. Plot the data points $$(x_i, y_i)$$, so we can see what it looks like. Clearly label your plot.

2. Using the linear least-squares methods discussed previously in class, find the coefficients of $$a_n$$ of the parabola $\hat{y}(x) := a_0 + a_1 x + a_2 x^2$ that best approximates the shape of the curve observed in our photo.

3. Make a new figure with 2 subplots. In the first subplot, show the data and your fit to the data, superimposed over each other. In the second subplot, plot the residuals $$[x_i, y_i-\hat{y}(x_i)]$$.

4. If your curve is a good fit, to all the data points, the residuals should be "independent", and should not show any patterns. Are your residuals random?

5. Calculate the total error $$E = \sum_{i} | y_i - \hat{y}(x_i) |$$ between your curve and the data.

6. In class, we predicted that the curve should be a catenary of the form $\tilde{y}(x) = b + c \cosh((x-d)/c).$ Try to find the $$(b,c,d)$$ that minimize the error between the catenary and the data. (Note: you can use your parabola coefficients and the observation that $$\cosh(x) \approx 1 + x^2/2$$ to find initial guesses for $$b$$, $$c$$, and $$d$$.

7. Compare the performance of Galileo's parabola hypothesis and Bernoulli's catenary hypothesis for explaining a hanging chain's shape. Discuss.

8. What about a general polynomial model of the form $\hat{y}(x) := a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 + a_5 x^5 + a_6 x^6 ?$

2. For a hanging chain, if we move the poles closer by half, how do the minimum and tension change?

3. Suppose that instead of having a uniform density, we have a chain whose density varies. Can we weight the chain so that it hangs in a circular arc rather than a catenary? If so, what would be the mass distrubition?

4. Suppose you building a simple rope bridge across a gorge in the Andes mountains. The paths are equal height on either side of the gorge. What length bridge should you build to minimize maximum force exherted on the bridge when it is empty and at rest?

5. Fermat's principle says that light rays travel the path of least-time between two points. While Fermat's principle is not quite right, it does explain common phenomena of reflection and refraction at once. Now, the refraction index of a material is define as the ratio of the speed of light in a vacuum to the speed of light in the material $$n = c / v$$.

Suppose we'd like to know the path of a light ray from the point $$(0,0)$$ to the point $$(1,0)$$ in a medium with refraction index $$n(x,y)$$.

$T = \int_{0}^1 \frac{n(x,y(x))}{c} \sqrt{1 + \left(\frac{dy}{dx}\right)^2} dx$

Use functional calculus to find a differential equation for the paths of light-rays between these two points.

6. (Bernoulli's Brachistochrome) Suppose a bead slides down a wire from an initial position $$(0,1)$$ to a final position $$(L,0)$$. Find a differential equation describing the shape of the wire which leads to the fastest fall of the bead between these two points.

7. F. Wan's derivation of the beam equation by functional differentiation. According to (Fung, 1965), (Sampson and the House of Dagon) Consider a vertical column of length $$L$$ compressed for a force $$P$$. If the potential energy of the Let $$x(y)$$ be the narrow displacement of this column from the perfect vertical at each height $$y$$, and $$\rho(y)$$ be the mass of the column per unit height. The potential energy of the beam is the sum of the bending energy and the gravitational potential energy minus the loss of energy from compression, or approximately $E[x] = \frac{1}{2} \int_{0}^L EI (x'')^2 - P (x')^2 dy$ Use functional differentiation to find an ordinary differential equation that will be satisfied when the potential energy is stationary (either being a minimum or a saddle point).

8. Suppose we have a large, thick piece of foam, that we compress allong a line using a narrow steel blade. Find an equation for the shape of surface of the foam around the compressed crease. You can assume all deformation is vertical and there is no lateral deformation.

9. Using an image of the Ead's bridge, extract a graph model of the bridge consisting of coordinates for nodes where steel beams meet and distances between connected nodes.

10. Although Galileo was incorrect about the shape of the hanging chain, his

11. Galileo's catenary: Show the equation for the forces on an elastic string under the pull of gravity, where the net force at each point is proportional to the second derivative in space leads to the 1-d PDE $m\frac{d^2u}{dt^2} = - \mu \frac{du}{dt} + \frac{d^2u}{dx^2} - mg.$