Catenary Between Two Points (Hanging Chain)

Question

This is one of several challenges left for the community by Calvin's Hobbies.

The curve that an idealised hanging rope or chain makes is a catenary.

^{Image by Bin im Garten, via Wikimedia Commons. Used under the CC-By-SA 3.0 licence.}

Write a program that will draw a catenary, as an image, in quadrant 1 of the plane given two points (x₁,y₁), (x₂,y₂), and the "rope length" L. L will be greater than the distance between the two points.

You must also draw axes on the left and bottom sides of the image (400x400 px min) for scale. Only draw the quadrant from x and y in range 0 to 100. (You may assume the points are in range.)

Dots, or something similiar, should be drawn at the (x₁,y₁), (x₂,y₂) endpoints to distinguish them. The curve should only be drawn in the space between these points.

Answer 1

Python + NumPy + Matplotlib, 1131

Just to get us started, here's an attempt that uses no knowledge of calculus or physics other than the fact that the catenary minimizes the energy of a chain. Hey, my algorithm may not be efficient, but at least it's not implemented efficiently either!

import math
import random
import numpy as np
import matplotlib.pyplot as plt
length, x0, y0, x1, y1 = input(), input(), input(), input(), input()
chain = np.array([[x0] + [length / 1000.]*1000, [y0] + [0.] * 1000])
def rotate(angle, x, y):
 return x * math.cos(angle) + y * math.sin(angle), -x * math.sin(angle) + y  * math.cos(angle)
def eval(chain, x1, y1):
 mysum = chain.cumsum(1)
 springpotential = 1000 * ((mysum[0][-1] - x1) ** 2 + (mysum[1][-1] - y1)  ** 2)
 potential = mysum.cumsum(1)[1][-1]
 return springpotential + potential
def jiggle(chain, x1, y1):
 for _ in xrange(100000):
  pre = eval(chain, x1, y1)
  angle = random.random() * 2 * math.pi
  index = random.randint(1,1000)
  chain[0][index], chain[1][index] = rotate(angle, chain[0][index], chain[1][index])
  if( pre < eval(chain, x1, y1)):
   chain[0][index], chain[1][index] = rotate(-angle, chain[0][index], chain[1][index])
jiggle(chain, x1, y1)
sum = chain.cumsum(1)
x1 = 2 * x1 - sum[0][-1]
y1 = 2 * y1 - sum[1][-1]
jiggle(chain, x1, y1)
sum = chain.cumsum(1)
plt.plot(sum[0][1:], sum[1][1:])
plt.show()

Answer 2

BBC Basic, 300 ASCII characters, tokenised filesize 260

  INPUTr,s,u,v,l:r*=8s*=8u*=8v*=8l*=8z=0REPEATz+=1E-3UNTILFNs(z)/z>=SQR(l^2-(v-s)^2)/(u-r)a=(u-r)/2/z
  p=(r+u-a*LN((l+v-s)/(l-v+s)))/2q=(v+s-l*FNc(z)/FNs(z))/2MOVE800,0DRAW0,0DRAW0,800CIRCLEu,v,8CIRCLEr,s,8FORx=r TOu
    DRAW x,a*FNc((x-p)/a)+q
  NEXT
  DEFFNs(t)=(EXP(t)-EXP(-t))/2
  DEFFNc(t)=(EXP(t)+EXP(-t))/2

Emulator at http://www.bbcbasic.co.uk/bbcwin/bbcwin.html

This has obviously been solved before, so the first thing I did was look what others have done.

The equation of a catenary centred at the origin is simply y=a*cosh(x/a). It becomes slightly more complicated if it is not centred at the origin.

Various sources say that if the length and endpoints are known the value for a must be determined numerically. There is an unspecified parameter h in the wikipedia article. So I found another site and basically followed the method here: http://www.math.niu.edu/~rusin/known-math/99_incoming/catenary

BBC Basic does not have sinh and cosh built in, so I defined two functions at the end of the program to calculate them using EXP

coordinates for lefthand point must be supplied before righthand point, OP confirmed this is OK. Length is given last. Values can be separated by commas or newlines.

Ungolfed code

  INPUT r,s,u,v,l

  REM convert input in range 0-100 to graphic coordinates in range 0-800 
  r*=8 s*=8 u*=8 v*=8 l*=8

  REM solve for z numerically
  z=0
  REPEAT
    z+=1E-3
  UNTIL FNs(z)/z>=SQR(l^2-(v-s)^2)/(u-r)

  REM calculate the curve parameters
  a=(u-r)/2/z
  p=(r+u-a*LN((l+v-s)/(l-v+s)))/2
  q=(v+s-l*FNc(z)/FNs(z))/2

  REM draw axes, 800 graphics units long = 400 pixels long (2 graphics units per pixel)
  MOVE 800,0
  DRAW 0,0
  DRAW 0,800

  REM draw markers at end and beginning of curve (beginning last, so that cursor is in right place for next step)
  CIRCLE u,v,8
  CIRCLE r,s,8

  REM draw curve from beginning to end
  FORx=r TOu
    DRAW x,a*FNc((x-p)/a)+q
  NEXT

  REM definitions of sinh and cosh
  DEF FNs(t)=(EXP(t)-EXP(-t))/2
  DEF FNc(t)=(EXP(t)+EXP(-t))/2

Answer 3

Python 2.7 + matplotlib, 424

Run as

python thisscript.py [length] [x0] [y0] [x1] [y1]

If I can assume that x0 is always smaller than x1 character count reduces to 398

from numpy import *
from pylab import *
from scipy.optimize import *
import sys
c=cosh
l,p,q,u,w=map(float,sys.argv[1:])
if p>u:
 p,q,u,w=u,w,p,q
h=u-p
v=w-q
a=brentq(lambda a:(2.*h/a*sinh(0.5*a))**2-l**2-v**2,1e-20,600)
b=brentq(lambda b:c(a*(1.-b))-c(a*b)-a*v/h,-600/a,600/a)
r=linspace(p,u,100)
plot([p,u],[q,w],'ro')
plot(r,h/a*c(((r-p)/h-b)*a)-h/a*c(a*b)+q,'k-')
gca().set_xlim((0,100))
gca().set_ylim((0,100))
show()

The magic number 600 you see appearing in some places is due to the fact that cosh(x) and sinh(x) start overflowing around x=710 (so 600 to keep some margin)

Basically I solve the problem in the frame where the catenary goes trough (0,0) and (x1-x0,(y1-y0)/(x1-x0)) and then remap to the original frame. This improves the numerical stability a lot.