[RESOLVED] catenary problem

[boops boops]

I'm trying to solve the classic problem of a (weightless) rope of length s hanging between any two points P1 and P2. Wikipedia gives the following guidance:



I can use the last formula to derive the value of a. My present method is to loop with trial values of a starting from 1 at intervals (e.g. 0.25) until the difference between the two sides of the equation passes through a minimum. No doubt there are more efficient ways to do the numerical analysis but it seems to work.

My problem is that, having derived the value of a, I cannot see how to relate the resulting catenary y = a * cosh (x/a) to the original locations P1 and P2. The Wiki formula depends only on the x difference (h) and y difference (v), not on the actual X and Y values. The instruction to "translate the axes" doesn't help me because I can't actually do that until I know the location of the vertex.

So can someone suggest how to "translate the axes" or an alternative approach?

BB

[jemidiah]

You want to solve for x_1, which will tell you how much things were translated horizontally; after that it's easy to figure out how much things were translated vertically. One way to solve for x_1 is to rewrite v and s in terms of exponentials. Let P = e^(x_1 / a), Q = e^(x_2 / a). Then...

v = a(e^(x_2 / a) + e^(-x_2 / a))/2 - a(e^(x_1 / a) + e^(-x_1 / a))/2
= a (Q + 1/Q - P - 1/P)/2
s = a(e^(x_2 / a) - e^(-x_2 / a))/2 - a(e^(x_1 / a) - e^(-x_1 / a))/2
= a (Q - 1/Q - P + 1/P)/2

=> (s+v)/a = Q - P
=> (s-v)/a = 1/P - 1/Q = (Q-P)/(QP) = (s+v)/a * 1/(QP)
=> QP = (s+v)/(s-v)
=> e^((x_2 + x_1) / a) = (s+v)/(s-v)
=> (x_2 + x_1) / a = ln( (s+v) / (s-v) )
=> x_2 + x_1 = a ln( (s+v) / (s-v) ) == A

Since x_2 - x_1 = h, we now have
x_1 = (A - h) / 2
= (a ln( (s+v) / (s-v) ) - h) / 2

You should double-check my derivation. There are probably simpler methods.

[boops boops]

Many thanks. Your formula for x_1 looks exactly what I need. I went through your derivation and I think I can just about understand it. I don't know if that counts as double-checking, but if I can implement it in code and conjure up a convincing-looking catenary on the screen it will be validation enough. I'll report back and resolve the thread if and when I get that far.

BB

[jemidiah]

I double-checked it myself and found I had forgotten the 1/2's in the exponential version of the hyperbolic trig functions. They canceled later though. I ran it through two edge cases and the results are at least sane: first, if everything is scaled by some factor, A is also scaled by that factor since it's proportional to a, which seems perfectly reasonable since the transcendental equation implies that a scales exactly by scale factors. Second, if v=0, the logarithm gives 0 since the inside becomes s/s=1, as it should since x_1 and x_2 will be on opposite sides of the origin. The real test cases are more involved and I won't take the time to try them myself, but in any case I've checked my formula enough that I'm confident in it. Mostly I was just sleepy when I did the derivation originally .

[boops boops]

I double-checked it myself and found I had forgotten the 1/2's in the exponential version of the hyperbolic trig functions. They canceled later though. I ran it through two edge cases and the results are at least sane: first, if everything is scaled by some factor, A is also scaled by that factor since it's proportional to a, which seems perfectly reasonable since the transcendental equation implies that a scales exactly by scale factors. Second, if v=0, the logarithm gives 0 since the inside becomes s/s=1, as it should since x_1 and x_2 will be on opposite sides of the origin. The real test cases are more involved and I won't take the time to try them myself, but in any case I've checked my formula enough that I'm confident in it. Mostly I was just sleepy when I did the derivation originally .

I had some difficulty cooking up a catenary with the X1 formula. It turned out, after much puzzling, that it was yielding the negative of x_2 instead of x_1, but I can't see why that happens from the derivation. Could I be making some different assumptions? Of course, it was easy to fix once I realized the problem: X_1 = - h - X_1.

For anyone interested, here's a VB.Net (Framework 4) demo of my catenary class: CatenaryDemo.zip You can drag the end points to any position, and change the length of the suspended "rope". Once again, many thanks to Jemediah.

BB

[jemidiah]

The assumptions are indeed slightly different, thanks to Windows' annoying inverted y-axis. If in Catenary_X1, pLeft.Y is 50 and pRight.Y is 100, pLeft is in fact above pRight by 50, so that v (== right y-value minus left y-value, assuming standard y-axis orientation) should be -50 whereas the line "Dim v As Double = pRight.Y - pLeft.Y" makes v=50 in this case. This flip has the effect of introducing a negative sign to each occurrence of v. The easiest fix is perhaps to use my original formula (which incidentally I am now entirely confident in, modulo obvious edge cases) while using "Dim v As Double = pLeft.Y - pRight.Y" in this particular routine. This is equivalent to your changes, though it's a little more direct. This wasn't a problem elsewhere mostly because v^2 was the important quantity and not v itself, so the sign change got washed out.

By the way, your changes work because...
(A(-v) - h)/2
= (a ln( (s+(-v)) / (s-(-v)) ) - h) / 2
= (a ln( (s-v) / (s+v) ) - h)/2
= (-a ln( (s+v) / (s-v) ) - h)/2
= -(a ln( (s+v) / (s-v) ) + h)/2
= -(A(v) + h)/2

P.S. Nice-looking hanging rope!

[boops boops]

The assumptions are indeed slightly different, thanks to Windows' annoying inverted y-axis. If in Catenary_X1, pLeft.Y is 50 and pRight.Y is 100, pLeft is in fact above pRight by 50, so that v (== right y-value minus left y-value, assuming standard y-axis orientation) should be -50 whereas the line "Dim v As Double = pRight.Y - pLeft.Y" makes v=50 in this case. This flip has the effect of introducing a negative sign to each occurrence of v. The easiest fix is perhaps to use my original formula (which incidentally I am now entirely confident in, modulo obvious edge cases) while using "Dim v As Double = pLeft.Y - pRight.Y" in this particular routine. This is equivalent to your changes, though it's a little more direct. This wasn't a problem elsewhere mostly because v^2 was the important quantity and not v itself, so the sign change got washed out.

By the way, your changes work because...
(A(-v) - h)/2
= (a ln( (s+(-v)) / (s-(-v)) ) - h) / 2
= (a ln( (s-v) / (s+v) ) - h)/2
= (-a ln( (s+v) / (s-v) ) - h)/2
= -(a ln( (s+v) / (s-v) ) + h)/2
= -(A(v) + h)/2

P.S. Nice-looking hanging rope!

Ah, I suspected the y inversion might have something to do with it but I didn't think it through. I'm thinking of posting a version of the class to the VB.Net codebank since I have not found anything similar on the web. Assuming you don't mind, it will contain your recommended X1 derivation with acknowledgements.

[jemidiah]

Assuming you don't mind, it will contain your recommended X1 derivation with acknowledgements.

Cool, go for it. I don't care about acknowledgment in this particular case since the derivation is (relatively speaking) simple--no insight (the hard part of doing math) required, just some fiddling with algebra. Thank you though for your thoroughness.

Cheers.

[loocas_duber]

Hi there, genltemen,

sorry to jump in, but this is probably the only thread I found that is really close to solving my issue.

I hope you're still around the forums, since I have trouble transforming your discussion into my code.

I still don't understand how to plot a curve (or several points laying on the curve) in 3D (or even 2D) from what you described here.

I get the formulas, but I don't get how to iterate over a set of points, or rather steps, at which I'd like to place these points along the curve.

Similarly I was dealing with bezier curves and easy enough, I just wrote an algo going from 0 to 1 and each step I placed a point laying on the bezier curve. After connecting the "dots" at small enough steps I got my bezier curve all nice and shiny in my 3D scene.

Now with this, all I'm getting from the calculations are single value outputs. Unfortunately I cannot understand VB code that much so I couldn't decipher boops' program I can program in Python, though, if it makes any difference

Anyways, all I need to do is to plot a curve (or a set of points laying on the curve) that is a freely hanging rope between two points in space.

The idea is that people specify the length of the rope/chain and pick the two points and the program will plot the catenaries between those two points based on the length of the curve. Obviously the length won't be allowed to be less than the distance between the two points.

So, my question is, how do I tranlsate the algorithms mentioned here into a loop that goes for as long as there are many steps (number of points laying on the curve/resolution of the curve) and get the X and Y coordinates of each of the point on the catenary?

Thank you very much in advance, best

- Lukas

[boops boops]

Hi Lukas,

Welcome to VBForums. This Math Forum can be a bit slow at times but it's mainly because people don't post questions very often. I think you'll find everyone is still here.

All the math(s) for the catenary is in the Catenary.vb class. You input the end points of the rope, the length of the rope and the point spacing you require, and it gives you back an array of points (CatenaryPoints) which represents the hanging rope. So it sounds like you need to concentrate on converting that code into Python or whatever language you prefer. I chose the variable names and commented the code in relation to the Wikipedia method in post #1 of this thread. Maybe that will help you understand it.

Out of curiosity, I tried a couple of online converters for vb.net to Python but none of them worked, even for the simplest snippets of code. I've no idea what was going wrong there. Maybe there is someone else on this forum who knows enough of both languages to give you some specific help with converting.

There is not much point in trying to convert the code in Form1.vb, which deals with painting the catenary points on the screen (size and colour etc.) and with mouse and trackbar interaction. I assume you can do that in your own Pythonesque way.

BB

[loocas_duber]

Hi Lukas,

Welcome to VBForums. This Math Forum can be a bit slow at times but it's mainly because people don't post questions very often. I think you'll find everyone is still here.

All the math(s) for the catenary is in the Catenary.vb class. You input the end points of the rope, the length of the rope and the point spacing you require, and it gives you back an array of points (CatenaryPoints) which represents the hanging rope. So it sounds like you need to concentrate on converting that code into Python or whatever language you prefer. I chose the variable names and commented the code in relation to the Wikipedia method in post #1 of this thread. Maybe that will help you understand it.

Out of curiosity, I tried a couple of online converters for vb.net to Python but none of them worked, even for the simplest snippets of code. I've no idea what was going wrong there. Maybe there is someone else on this forum who knows enough of both languages to give you some specific help with converting.

There is not much point in trying to convert the code in Form1.vb, which deals with painting the catenary points on the screen (size and colour etc.) and with mouse and trackbar interaction. I assume you can do that in your own Pythonesque way.

BB

Hey boops,

thanks for letting me know. Yeah, I went through your code and, unfortunately, I got frustrated after a while of trying to decipher the variables and the general logic of the code.

The thing is, I'm not going to end up writing the Catenary plotter in Python, I just thought that it's such a popular language that we can discuss the technicalities directly in it, if need be. But I'm rather more interested in the actual algos in obtaining the points laying on the catenary.

The problem, for me, is that I really don't know where to start. I have a bunch of math equtions that I kind of understand. But how do I translate them to my specific programming language? I need to create a loop of a specific number of iterations and for each loop I need to find X and Y coordinates of a point that lays on the catenary that is hung between P1 and P2 in space (say 2D space, as 3D is just another transformation on top of the 2D one).

Thanks again, cheers...

[boops boops]

If you have specific questions about individual functions in the catenary.vb file, I could probably give you some further explanation. Unfortuantely I can't help you about porting the code to Python. Hang on, though, maybe someone else can.

BB

[loocas_duber]

If you have specific questions about individual functions in the catenary.vb file, I could probably give you some further explanation. Unfortuantely I can't help you about porting the code to Python. Hang on, though, maybe someone else can.

BB

Thanks, no problem really.

I might have figured it out, though. I'll have to test some more to confirm my "suspicion"

Thank you again, cheers.

[jemidiah]

@loocas_duber: It's probably best for you to get very basic functionality working first, and then add on other stuff. In particular, you should get something that can just plot y = cosh(x) in 3D. After that you can make things fancier by allowing the user to specify a pair of points and draw the catenary through them. The above algebra and discussion gives the vast majority of the ideas needed for that, though.

[loocas_duber]

@loocas_duber: It's probably best for you to get very basic functionality working first, and then add on other stuff. In particular, you should get something that can just plot y = cosh(x) in 3D. After that you can make things fancier by allowing the user to specify a pair of points and draw the catenary through them. The above algebra and discussion gives the vast majority of the ideas needed for that, though.

Hey, jemidiah, thanks for jumping in!

So, I've managed, from a few various discussions on the internet as well as with the Wikipedia help, to write a function in my 3D app, that plots a catenary using point helper objects (simple crosses in the 3D viewport). The problem I'm having now is the accuracy of the curve. It seems there's something wrong.

First things first:

I define a Line Length (the catenary rope length), s, then I derive the horizontal distance from the two points (P1 and P2) the catenary should be hung on, d = P2.x - P1.x, assuming P2 is farther on the X axis than P1. Then, I define the vertical distance between the points, h = P2.z - P1.z (assuming the P2 is higher on the Z axis, than P1. Also note that my 3D software considers the Z axis as the "up" axis). Then I have to find the parameter a. I wrote a function that takes the a, d, s, h, epsilon (for error margin checking) and iterationSteps and computationally try to find the a, based on the Newton-Raphson method:

i = 0
tempResult = a
result = squareRoot(s^2 - h^2)

while i <= iterationSteps do
(
fx = 2 * tempResult * sinh(d/(2 * tempResult))
dFx = 2 * sinh * (d/(2 * tempResult)) - (d * cosh(d/(2 * tempResult))/tempResult) #basically just a derivative of fx

if abs(fx) - result <= epsilon do exit #this is a test to see if the end result of the calculation is close enough to the actual result of the function based on an error margin epsilon

tempResult -= fx/dfx

i += 1
)

This function should return an approximately very close value of a for the catenary function.

Now, for the actual calculation of the Catenary curve between ponts P1 and P2, which, as I read, has to be calculated with a specific offset in order to fit the curve between the points. Here's my approach:

x1 = (a * log((s + h)/(s-h)) - d)/2

for the actual X offset:

xOffset = x1 - P1.x
zOffset = P1.z - a * cosh(x1/a)

This is all I need in order to calculate the actual Catenary curve. All seems to be working fine. Here's the last piece of the puzzle, the actual plotting of points that lay on the Catenary curve:

for x = 0 to numberOfPoints do
(
pointXPosition = ((P2.x - P1.x)/numberOfPoints) * x
pointZPosition = a * (cosh(pointXPosition + xOffset)/a) + zOffset
)

So far so good. The function really does plot the catenary and for the test example I was drawing on, it works flawlessly. Here are the parameters for my text example: P1=(0,10), P2=(20,15) d=20, h=5, s=20.7, the equation 2a*sinh(d/2a) = sqrt(s^2-h^2) yields a value for parameter a = 61.9388. In the end, the curve looks just like the guy's from the forums I drew on:



But as soon as I change any of the input parameters (position of points P1 and P2 or the lenght of the catenary rope) the result is just wrong!

Here, the P2 was moved slightly up the Z axis and the length was changed to 30.7, this is what I got:


Here, the P2 was moved closer to P1 and the rope length was changed to 60.0, this is the result:


Is there anything I've overlooked? Or what could be the cause for this errornous behavior?

Thanks a lot in advance, any help will be much appretiated!

[jemidiah]

Your Newton-Raphson code seems to have a lot of errors. Most of your code is actually (though incorrectly) computing a root for 2a sinh(d/(2a)), but "if abs(fx) - result <= epsilon" is masking the errors at least a little--presumably "result" is getting larger than abs(fx) eventually (you are finding a root of fx after all), making the difference negative and trivially less than epsilon, allowing you to at least exist the loop. Your derivative is also incorrect as written; I don't know why "sinh" appears by itself.

Hopefully the following corrects all the errors:

i = 0
tempResult = a
result = squareRoot(s^2 - h^2)

while i <= iterationSteps do
(
fx = 2 * tempResult * sinh(d/(2 * tempResult)) - result
dFx = 2 * sinh(d/(2 * tempResult)) - (d * cosh(d/(2 * tempResult))/tempResult) #basically just a derivative of fx

if abs(fx) <= epsilon do exit #this is a test to see if the end result of the calculation is close enough to the actual result of the function based on an error margin epsilon

tempResult -= fx/dfx

i += 1
)

You should make sure the output of this method, whatever it is, actually does satisfy 2a sinh(d/(2a)) ~= sqrt(s^2 - h^2) for some test inputs. I don't have time to look through the rest at the moment, though hopefully these changes will help. There are also a number of stylistic problems with your code, but the same is true of everyone when they start coding. (Sorry if this post seems negative, I'm a bit of a perfectionist.)

[loocas_duber]

Thank you very much, jemidiah! I appretiate it!

I'll test your code later today and post my results.

As for stylistic issues with my code, I'd very much appretiate, if you could point out concretely so I can learn from the mistakes. I, obviously, want to write clean and functional code, so any type of correction or pointing out mistakes from your or anyone else's side will be much appretiated.

Thank you again, I'll let you know how it went.

I don't know why "sinh" appears by itself.

Sorry, that was a typo, the multiplier wasn't supposed to be there. I took the derivative of the function from Wolfram as I wouldn't be able to derive it by hand myself so I assume it is correct.

[loocas_duber]

So, jemediah, unfortunately subtracting the sqrt(s^2-b^2) from the result of fx first (your suggested modification), didn't make any difference

[jemidiah]

Here's a Python translation of my previous code. It gives correct output for a(10, 5, 3). Just to be clear, I modified more than just the fx line. Does your version calculate a correctly while still producing incorrect output later?

def a(s, h, d):
 i=0
 tempResult = 1
 result = math.sqrt(s**2 - h**2)
 while i <= 500:
  fx = 2 * tempResult * math.sinh(d/(2 * tempResult)) - result
  dFx = 2 * math.sinh(d/(2 * tempResult)) - (d * math.cosh(d/(2 * tempResult)) / tempResult)
  if abs(fx) <= 0.0001: break
  tempResult -= fx/dFx
  i += 1
 return tempResult, fx, result

Edit: As for stylistic issues, some that come to mind immediately...
(Note: You didn't post your full code, and it seems you didn't just copy and paste it, so some of this may not be relevant.)

* Bad variable names--tempResult is annoyingly long to type yet not terribly informative; "a" would probably be better. "result" is probably more intuitively named "target" or simiar. "dFx" capitalizes the F for no apparent reason, and the x isn't helpful. "df" and "f" are probably better. "epsilon" is probably alright, though "errorMargin" is more informative and fits with the iterationSteps naming convention without being much longer.

* Bad comments: what you did comment wasn't very efficient/helpful, and you didn't write the important things. You're using Newton's method for root finding, which should probably be at the top of the routine; let someone Google it for an explanation. "basically just a derivative of fx" is incorrect on two counts--it's not just "basically" a derivative, it is the derivative. The long comment on the "if" line is just plain unhelpful--anyone who can read code will understand what that line is doing (well, if it were written correctly). In writing comments, you want to give design details that can't be instantly gleaned from just reading the code.

* Not very clean: your Newton's Method code could be made cleaner in a few ways. Switching from "while" to "for" would get rid of the "i += 1" at the end, the "i=0" at the beginning, and be more intuitive. You repeatedly use (d/(2 * tempResult)) in your formulas, so it might be nice to put this in a variable by itself, just for readability. (A compiler should optimize out the inefficiency of recomputing it anyway, so that shouldn't be an issue.)

* Potentially bad design: there are a zillion implementations of Newton's Method or more powerful numerical solvers out there. Depending on your language, it might be best to just import a numerical library and have it do the work. The main benefit of this general approach IMO is that someone else has done the debugging work for you. Also, depending on the language, there are more intuitive implementations of Newton's Method: you could perhaps create a general routine that takes in a pair of functions as an argument, the target function and the derivative, then runs the algorithm on it.

* Whitespace use: indentation is very helpful, always indent properly, it's worth the time. Your long comment spans multiple lines and should instead be broken up. A lot of programmers use ~80 characters per line, though this varies. There's also a lot of blank lines in your code which may or may not belong there--I'd have to see how the rest of the code looks to have a solid opinion.

[loocas_duber]

Thank you, jemidiah, I'm really grateful for your input. I'll have to investigate what's causing the trouble. The Newton method function seems to be producing "OK" results, i.e. close enough, but still not as close as when writtein in Pyhon, or computed manually, which is strange. One issue in my particular language of choice is the fact that all the trigonometric functions assume as input DEGREES, while in Python, the same functions assume RADIANS, which is a bit problematic as I always have to keep an eye on what is what and convert accordingly. That might be one of the issues (though I've tried numerous variations, with the same errornous results).

[jemidiah]

For the trig function issue, you can always define new functions like "cosr(x) = cos(radToDeg(x))". That could help take off some of the "cognitive load" (roughly, the things you have to keep track of). I haven't heard of hyperbolic trig functions taking degrees but I suppose it happens....

What do you mean by "close enough"? For instance, what value of "a" does your Newton's Method give for s=10, h=5, d=3?

[loocas_duber]

I'm getting a = 57.2958 for the values you provided.

[jemidiah]

The value my routine gives with s=10, h=5, d=3 is a = -0.5396923312383373. Here sqrt(s^2 - h^2) = 8.660254037844387 while with my a, 2*a*math.sinh(d/(2*a)) = 8.660254751809196, which is very close. With your value of a, I get 2*a*math.sinh(d/(2*a)) = 3.000342706096323, which is very far from correct.

[loocas_duber]

Now that is a problem. But where do I start if the algo is correct (assuming it is, as I simply translated wikipedia's Newton's method to my scripting language).

EDIT: Got it! It was a problem of the degrees to radians and vice versa conversions.

Thanks a million, jemediah! I appretiate your help.

By the way, what would be your approach for solving the same curve, but the second point would be behind the first one? I'm going to try to work around this issue by checking the second point first, if it's to the left (behind) the first point, I'm going to make it the first point instead. Or, is there a more elegant solution to this issue?

[jemidiah]

Your proposed solution sounds great. It reuses your existing work and only adds a small amount of complexity to the logic. There might be a more efficient way, but I won't take the time to look for one.

[loocas_duber]

Sure thing, I don't want to bother you with this anymore.

Just so you know, by the way, I found even better solution, in my opinion, which will allow me to plot the catenaries in three dimensions: I'll construct a local world transform originating at the first points position and aligned to the second point with the global Up axis as the up vector. That way all my calculations will stay the same, only at the end all the points laying on the catenary curve will get transformed back to the original world transform according to the local world transform relative transform to the global world. If that makes any sense This way I won't have to deal with what point is behind or in front as the point number 2 will always be in front of the point no. 1 and also, I'll be able to transform the catenary curve any way I like in the 3D world. A bit more complicated, but a lot more useful.

Anyways, thank you again, you've been most helpful! Best.

[jemidiah]

That sounds like a fine solution. You can even confuse people and change the direction of gravity

Glad to have been helpful.