Modelling Question
Further readings
Bukowski, 2008 on Christian Huygens study of the hanging chain.
Clifford Truesdell, The Rational Mechanics of Flexible or Elastic Bodies, 1638-1788, 1960, also on the hanging chain.
A First Course in the Calculus of Variations by Mark Kot, 2014,has a beautiful collection of case studies.
Variational Principles of Mechanics (free) by Cornelius Lanczos, 1949 is a very readable classic text.
Introduction to the Mathematical Theory of Control by Bressan and Piccoli, 2007, is a concise modern formal treatment.
Lecture Notes on the Calculus of Variations (free) by Oskar Bolza,1904
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.
There are several different approaches for trying to come up with an equation for the shape of a hanging chain.
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{k}{\rho g} \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} \right)^2 = 1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2 \end{gather*}\] 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.
Despite being a nonlinear equation, we are fortunate that this equation may be integrated directly. We can make our calculation a little more compact by introducing a new constant \(H = \rho g/k\). \[\begin{gather} \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} = H \sqrt{ 1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2 } \\ u = y', \quad u' = y'' \\ \int \frac{ \mathrm{d} u }{ \sqrt{ 1 + u^2 } } = H \int \mathrm{d} x \\ u = \sinh \theta, \quad \mathrm{d} u = \cosh \theta \mathrm{d} \theta, \quad \cosh^2\theta - \sinh^2 \theta = 1 \\ y' = u = \sinh \theta = \sinh( g \rho (x + C) ) \\ y(x) = \frac{1}{H} \cosh( H (x + C) ) + K \end{gather}\] We can also verify the solution using sympy
from sympy import *
from sympy.abc import x,y,m,a,c,w
P = Integral( y(x) * sqrt( y(x).diff(x) ** 2 + 1 ), (x,0,w))
L = Integral( sqrt( y(x).diff(x) ** 2 + 1 ), (x,0,w))
A = P.args[0]-m * L.args[0]
F = (A.diff(y(x)) - A.diff(y(x).diff(x)).diff(x)).factor()
F.subs(y(x), m+cosh(c * (a+x))/c).doit().factor().simplify()
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\) and \(K = -\cosh(Hw/2)/H\). \[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 satisfy two 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\).
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.
(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.
A simple implementation, using a sequence of nodes on the string.
from numpy import array, arange, linspace
from scipy import sqrt
from matplotlib.pyplot import figure,show,xlim,ylim,title
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
dl = 0.2/N # resting distance
rho = 1./N # node mass
k = 4.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 ( i -1, i + 1 ):
v = x[j] - x[i], y[j] - y[i]
dist = sqrt(v[0]**2 + v[1]**2)
v = v[0]/dist, v[1]/dist
D[2*N+i] += v[0]*k*(dist-dl)
D[3*N+i] += v[1]*k*(dist-dl)
# 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 = 8
dt = 0.01
time=arange(start,end,dt)
X = array([0.]*(4*N))
X[0:N] = linspace(0,1,N)
Y = integrate.odeint(vfield,X,time)
def animate(i):
x = animate.U[i,0:N]
y = animate.U[i,N:2*N]
animate.line.set_data(x,y)
animate.title.set_text('Time t=%03.2f'%animate.time[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.U = Y
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()
Consider usinga a potential-formulation for force-distance calculations.
Use a graph-model, which provides a more general solution and can accomodate various truss designs and node masses.
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()
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
Plot the data points \((x_i, y_i)\), so we can see what it looks like. Clearly label your plot.
Using the linear least-squares methods discussed previously in class, find the coefficients of \(a_n\) of the polynomial \[\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\] that best approximates the shape of the curve observed in our photo. List all 7.
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)]\).
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?
Calculate the total error \(E = \sum_{i} | y_i - \hat{y}(x_i) |\) between your curve and the data.
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. Discuss.
For a hanging chain, if we move the poles closer by half, how do the minimum and tension change?
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?
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?
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.
(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.
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).
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.
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.