Take the 2-minute tour ×
Physics Stack Exchange is a question and answer site for active researchers, academics and students of physics. It's 100% free, no registration required.

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?

share|improve this question
    
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

1 Answer 1

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];
}
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.