Catenary using N straight lines

Question

I would need to draw catenary between dots A and B, A(x1, y1, z1), B(x2, y2, z2) using N straight lines.

I tried with formula for catenary

z = a * Math.cosh(x / a)

where

a = H / (H = horizontal tension, = weight per unit)

and my code is

xCurrent = x1;
yCurrent = y1;
zCurrent = z1;

for (N number of lines) {
  x1 = xCurrent;
  x2 = xCurrent + x1x2TotalDistance/N;
  y1 = yCurrent;
  y2 = yCurrent + y1y2TotalDistance/N;

  H = 25000;
  omega = 1.4 * length;
  a = H / omega;

  z1 = a * Math.cosh(x1 / a);
  z2 = a * Math.cosh(x2 / a);
  zDifference = z1 - z2;
  zNew = zCurrent - zDifference;

  DrawLine(T1(x1, y1, zCurrent), T2(x2, y2, zNew))

  zCurrent = zNew;
  xCurrent += x1x2TotalDistance/N;
  yCurrent += y1y2TotalDistance/N;
}

But the problem is that last zNew is not equal to z2 from B dot. Difference is not that big but it is too big to be noticable.

And if I increase H for example 100 times, then last zNew is equal to z2 in acceptable way, but I don't have catenary shape anymore, it draws straight line.

Is there a way to have last zNew equal to z2 from B dot to at least 2 decimal points, and to keep shape of catenary?

Answer 1

Let's assume the two points are $\vec{p}_1 = (x_1 , y_1 , z_1)$ and $\vec{p}_2 = ( x_2 , y_2 , z_2 )$, with acceleration due to gravity towards negative $z$ axis. Then, we have $$\left\lbrace ~ \begin{aligned} x(t) &= (1-t)x_1 + t x_2 = x_1 + t (x_2 - x_1) \\ y(t) &= (1-t)y_1 + t y_2 = y_1 + t (y_2 - y_1) \\ z(t) &= z_0 + \displaystyle a \cosh\left( \frac{w (t - t_0)}{a} \right) \\ \end{aligned} \right.$$ where $w$ is the horizontal distance, $w = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)}$, $t_0$ is the relative location of the lowest point of the catenary (to be solved), $z_0$ is the $z$-location of the catenary (also to be solved), and $a$ is the scale factor: $$a = \frac{H}{m}$$ where $H$ is the horizontal tension, and $m$ is the total mass of the catenary.

To determine $z_0$ and $t_0$, we need to apply $$\left\lbrace ~ \begin{aligned} z(0) &= z_1 \\ z(1) &= z_2 \\ \end{aligned} \right.$$ Subtracting the latter from the former gives us $$\cosh\left(\frac{w (1 - t_0)}{a}\right) - \cosh\left(\frac{w t_0}{a}\right) - \frac{z_2 - z_1}{a} = 0$$ which (I believed!) we have to solve numerically for $t_0$.

In practice, we can use $$f(x) = \cosh\bigr(k (1 - x)\bigr) - \cosh\bigr(k x\bigr) - c$$ where $k = w/a \gt 0$ and $c = (z_2 - z_1)/a$. $f(x)$ is monotonically decreasing everywhere. So, we can start at $x = 1/2$, and extend the search right or left depending on the sign of $f(x)$. When the range of $x$ spanning zero is found, you can use the binary search method on that span to find $x$.

Here is a Python 3 example for finding $t_0 = x$ given $k$ and $c$:

from math import sinh, cosh

def find_t0(k, c, epsilon = 0.0000001):
    # Find xmin..xmax where f(xmin) > 0 and f(xmax) < 0
    if c < 0.0:
        step = 1.0
        xmin = 0.5
        xmax = xmin + step
        while cosh((1.0 - xmax)*k) - cosh(k*xmax) - c > 0.0:
            xmin = xmax
            xmax = xmax + step
            step = step * 2
    elif c > 0.0:
        step = 1.0
        xmax = 0.5
        xmin = xmax - step
        while cosh((1.0 - xmin)*k) - cosh(k*xmin) - c < 0.0:
            xmax = xmin
            xmin = xmin - step
            step = step * 2
    else:
        # Special case, c = 0, so f(0.5) = 0.
        return 0.5

    # Use a binary search to find x.
    while xmax - xmin > epsilon:
        x = (0.5 * xmax) + (0.5 * xmin)
        f = cosh((1.0 - x)*k) - cosh(k*x) - c
        if f > 0.0:
            xmin = x
        elif f < 0.0:
            xmax = x
        else:
            return x

    return (0.5 * xmax) + (0.5 * xmin)

Note that epsilon is the desired relative precision; i.e. 0.000001 gives you about six digits of precision.

Edited: Claude Leibovici showed a proper mathematical solution in their answer here. So, instead of the above, I recommend using

from math import sinh, cosh, log

def find_t0(k, c):
    if c == 0.0:
        return 0.5

    a = cosh(k)
    b = sinh(k)

    # This looks odd, but we need to calculate d = 1 + b - a
    # this way to avoid domain cancellation error.
    d = 1.0 - (a - b)

    r = sqrt(c*c + b*b - a*a + a + a - 1.0)
    # Since r > |c| for k > 0, (r - c)/d > 0 too
    return log((r - c) / d) / k

I have verified with a large number of random $k \gt 0$ and $c$ that the two implementations produce the same results (within the specified epsilon for the numerical search) – but this latter version is obviously much, much more efficient! (This even better version now without any conditionals, again thanks to Claude Leibovici.)

The only trick there is that for $k \ge 39$ or so, $\cosh(k) - 1 = \cosh(k) = \sinh(k)$ using double-precision floating-point numbers. This means that this operation is suspect to domain cancellation error. Fortunately, $d(k) = 1 + \sinh(k) - \cosh(k)$ is a simple function with $d(0) = 0$, $d(k) \lt 1$ for $k \in \mathbb{R}$, and we can avoid the domain cancellation error by subtracting the two hyperbolic functions first, then adding 1 to their difference. (In C, one would use e.g. d = 1.0 - (double)(a - b);, as the cast forces a - b to be evaluated first and at double precision and range.) This will be as accurate as is possible with floating point numbers for the entire range $k \gt 0$, $k \in \mathbb{R}$ we're interested in.

Finally, we can solve $z_0$ from $z(0) = z_1$ (or $z(1) = z_2$), i.e. $$z_0 = z_1 - a \cosh\left( \frac{- w t_0}{a} \right)$$

To draw the catenary using $N$ line segments $(x_i, y_i, z_i)$ to $(x_{i+1}, y_{i+1}, z_{i+1})$, for $i = 0 .. N-1$. Each point $(x_i, y_i, z_i)$ for $i = 0 .. N$ is $$\begin{aligned} t &= \frac{i}{N} \\ x_i &= (1-t)x_1 + t x_2 \\ y_i &= (1-t)y_1 + t y_2 \\ z_i &= z_0 + \displaystyle a \cosh\left( \frac{w (t - t_0)}{a} \right) \\ \end{aligned}$$

Answer 2

Starting from @Guest' answer $$f(x) = \cosh\bigr(k (1 - x)\bigr) - \cosh\bigr(k x\bigr) - c$$ use $$ \cosh\bigr(k (1 - x)\bigr)=\cosh (k) \cosh (k x)-\sinh (k) \sinh (k x)$$ Let $a=\cosh (k)$, $b=\sinh (k)$, $y=kx$ to face $$g(y)=(a-1) \cosh (y)-b \sinh (y)-c$$ Now, let $y=\log(t)$ and you need to solve $$(1+b-a)t^2+2 c t+(1-a-b) =0$$ for which $$\Delta=2(\cosh(k)-1)+c^2 \qquad > 0\qquad \forall k \quad \text{and} \quad \forall c$$

The product of the roots being $-2(\cosh(k)-1) <0$, pick the positive solution $t_+$ and then $y=\log(t_+)$ and finally $x=\frac{ \log(t_+)}k$ .

$$x=\frac 1 k \log \left(\frac{c-\sqrt{c^2+2 \cosh (k)-2}}{\cosh (k)-\sinh (k)-1}\right)$$