# Closed form for shape/tension of an elastic cable slung between two points

Given the 2D coordinates of two points, $a$ and $b$, between which an elastic cable of known length, $l$, mass per unit length, $m$, and the spring constant, $e$, is slung, I need to compute the shape of the cable, and also the horizontal tension, $t$, in the cable.

So far I have the equations for the x and y coordinates of the cable, parameterized by $p$, which is the distance along the unstretched cable: $$f_x(p) = \frac{t}{mg} \sinh^{-1}\left(\frac{mgp}{t}\right) + \frac{tp}{e} + c_x.\\ f_y(p) = \sqrt{\left(\frac{t}{mg}\right)^2 + p^2} + \frac{mgp^2}{2e} + c_y.\\$$ where $c_x$ and $c_y$ are constants, which leads me to the following triplet of simultaneous equations: $$f_x(q) - f_x(r) = a_x - b_x.\\ f_y(q) - f_y(r) = a_y - b_y.\\ |q-r| = l.$$ in three unknowns, $t$, $q$ and $r$ (given that the constants cancel), where $q$ and $r$ are the values of the parameter $p$ at the points $a$ and $b$ respectively.

How would you compute those unknowns, and can it be done in closed form?

-
Does $t$ depend on the position? Is it the tension at some extremity? – fffred Jul 2 '13 at 20:34
Related: physics.stackexchange.com/q/51485/2451 and links therein. – Qmechanic Jul 3 '13 at 0:13
fffred: No, the value $t$ is constant along the cable. It is the $x$ (horizontal) component of the tension at any point in the cable. – user664303 Jul 3 '13 at 8:54
Qmechanic: Thanks. I googled "catenary" before posting, and this lead to what I have already. – user664303 Jul 3 '13 at 8:59

I don't know about a closed form solution, but I use the Gauss-Newton optimization approach to solve the problem, in fewer than 20 iterations usually, using the C code below.

#include <cmath>

/* Automatically generated functions to compute residuals and Jacobian of residuals */
inline void compute_residuals(const double data[5], const double params[2], double res[2])
{
double t2, t3, t4, t5, t6, t7, t8, t9;
t2 = 1.0/data[3];
t3 = 1.0/params[0];
t4 = params[1]+data[2];
t5 = params[0]*params[0];
t6 = 1.0/(data[3]*data[3]);
t7 = t5*t6;
t8 = params[1]*params[1];
t9 = t4*t4;
res[0] = -data[0]+t2*params[0]*log(sqrt((t3*t3)*(t4*t4)*(data[3]*data[3])+1.0)+t3*t4*data[3])-t2*params[0]*log(sqrt((t3*t3)*(params[1]*params[1])*(data[3]*data[3])+1.0)+t3*params[1]*data[3])+t4*params[0]*data[4]-params[0]*params[1]*data[4];
res[1] = -data[1]-sqrt(t7+t8)+sqrt(t7+t9)-t8*data[3]*data[4]*(1.0/2.0)+t9*data[3]*data[4]*(1.0/2.0);
}
inline void compute_jacobian(const double data[5], const double params[2], double J[4])
{
double t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20, t21, t22, t23, t24;
t2 = params[1]+data[2];
t3 = 1.0/(params[0]*params[0]);
t4 = data[3]*data[3];
t5 = 1.0/data[3];
t6 = 1.0/params[0];
t7 = params[0]*params[0];
t8 = 1.0/(data[3]*data[3]);
t9 = t2*t2;
t10 = t7*t8;
t11 = params[1]*params[1];
t12 = t3*t4*t9;
t13 = t12+1.0;
t14 = 1.0/sqrt(t13);
t15 = t3*t4*t11;
t16 = t15+1.0;
t17 = 1.0/sqrt(t16);
t18 = t10+t11;
t19 = 1.0/sqrt(t18);
t20 = t9+t10;
t21 = 1.0/sqrt(t20);
t22 = params[1]*2.0;
t23 = data[2]*2.0;
t24 = t22+t23;
J[0] = -t2*t14+t17*params[1]+t5*params[0]*log(sqrt((t2*t2)*(t6*t6)*(data[3]*data[3])+1.0)+t2*t6*data[3])-t5*params[0]*log(sqrt((t6*t6)*(params[1]*params[1])*(data[3]*data[3])+1.0)+t6*params[1]*data[3])+t2*params[0]*data[4]-params[0]*params[1]*data[4];
J[1] = -t7*t8*t19+t7*t8*t21;
J[2] = t14-t17;
J[3] = t21*t24*(1.0/2.0)-t19*params[1]+t24*data[3]*data[4]*(1.0/2.0)-params[1]*data[3]*data[4];
}

/* Update the parameters with an iteration of Gauss-Newton, and return the score of the starting position */
inline double gauss_newton_iteration(const double data[5], double params[2])
{
/* Compute the residuals and Jacobian of resdiduals */
double score, res[2], J[4], H[3], det;
compute_residuals(data, params, res);
compute_jacobian(data, params, J);

/* Compute the score, inverse Hessian and Jacobian of the cost function */
score = res[0] * res[0] + res[1] * res[1];
H[2] = J[0] * J[0] + J[1] * J[1];
H[1] = -(J[0] * J[2] + J[1] * J[3]);
H[0] = J[2] * J[2] + J[3] * J[3];
det = 1.0 / (H[0] * H[2] - H[1] * H[1]);
H[0] *= det;
H[1] *= det;
H[2] *= det;
J[0] = J[0] * res[0] + J[1] * res[1];
J[1] = J[2] * res[0] + J[3] * res[1];

/* Compute the parameter update */
params[0] *= exp(-(H[0] * J[0] + H[1] * J[1]));
params[1] -= H[1] * J[0] + H[2] * J[1];
return score;
}

/* Given an initial Pa, compute T0 */
inline double compute_T0(const double data[5], const double Pa)
{
double T0 = Pa + data[2];
T0 = data[1] + 0.5 * data[3] * data[4] * (Pa * Pa - T0 * T0);
T0 = (T0 + data[2]) * (T0 - data[2]) * (data[2] - T0 + 2.0 * Pa) * 0.5 / (T0 + 1e-300);
T0 = sqrt(abs(T0));
return T0 * data[3] + 1e-15;
}

/* Given parameters and distance along the line (0 being start, 1 being end), compute a tension */
inline double compute_tension(const double data[5], const double params[2], double P)
{
P = params[1] + data[2] * P;
P = sqrt(P * P * data[3] * data[3] + params[0] * params[0]);
return P;
}

/* Given a 3D start point and end point, compute the tension at those two
* points of an elastic line slung between the points.
*
*IN:
*   start - 3D coordinates (XYZ in metres) of the start point of the line (gravity acts in the -Z direction).
*   end   - 3D coordinates (XYZ in metres) of the end point of the line (gravity acts in the -Z direction).
*   l     - scalar length of line (unstretched) between the start and end points.
*   m     - mass per unit length (in Kg/m) of the line (unstretched).
*   e     - coefficient of elasticity of the line.
*   Pa    - estimate of position along the line (in m, unstretched) of the start point, zero being the horizontal point of the catenary.
*
*OUT:
*   T     - line tensions (in N) at start, end, and the point of zero vertical tension.
*   return value - computed value of Pa.
*/
double elastic_catenary(const double start[3], const double end[3], const double l, const double m, const double e, const double Pa, double T[3])
{
int iter;
double last_score, score = 1.0e300;
double params[2];
double data[5] = {0, end[2]-start[2], l, m * 9.81, 1.0 / e};
data[0] = sqrt((end[0] - start[0]) * (end[0] - start[0]) + (end[1] - start[1]) * (end[1] - start[1]));

/* Initialize parameters */
params[1] = Pa;
params[0] = compute_T0(data, Pa);

/* Do the optimization */
for (iter = 0; iter < 100; ++iter) {
last_score = score;
score = gauss_newton_iteration(data, params);
if (score <= last_score && (last_score - score) <= (last_score * 1e-10))
break; /* Converged */
}

/* Compute the tensions */
T[0] = compute_tension(data, params, 0); /* Tension at line start */
T[1] = compute_tension(data, params, 1); /* Tension at line end */
T[2] = params[0];                        /* Tension at the horizontal point  */

/* Return computed value of Pa - useful for initializing future optimizations */
return params[1];
}

-