Numerical approximation method needed for solving hanging rope equation system

[SevenSeas]

Sorry, this is not really algebra but I didn't find a more suitable forum. What I'm in need of is a suitable numerical approximation method for solving the equation system below, a numerical approximation method that I then can implement in javascript for myself.

As you may know from another thread I posted in Calculus recently, but which I found the answer to myself, I'm writing a javascript game about hanging in threads.

Now, when an idealised thread of length L isn't stretched, it will be hanging in an arc between two points of heights H (higher height) and K (lower height).

I found almost all the information I needed to implement such a thread model on this page: http://www.had2know.com/academics/catenary-equation-shape-hanging-chain.html

It tells the following:

The equation for the hanging thread is as follows:

y(x) = a + (1/b)cosh(b(x-c))

where a, b, and c are constants satisfying the following equations:

1. H = a + (e^bc + e^-bc)/(2b)
2. K = a + (e^b(D-c) + e^-b(D-c))/(2b)
3. L = (e^b(D-c) - e^-b(D-c) + e^bc - e^-bc)/(2b)

and as told, we know the values of H, K and L.

But unfortunately, this equation system is not solvable algebraically, but must be solved numerically. And I know very little about methods for numerical approximation. And I won't be able to plug Matlab or something similar into my code, but I will have to write in computer code some (hopefully simple) approximation function myself, suitable for solving this system with three unknown variables.

There was even a calculator for just this on the page I linked, but as I understand it the code for that calculator seems to reside on the server, so I can't just "be inspired" by that calculator's code...

I checked the Wikipedia page about numerical approximation but it was hard to sort out all the information.

So, do you have a suggestion of which method to use here? Thank you in advance!

 

[SevenSeas]

Oh, and the constant D is the distance between the two poles from which the rope is hanging, forgot to say that.. :cool:

I'm checking out Newton's method generalised for equation systems of n variables now, trying to calculate all the partial derivatives of these functions above, for the Jacobi matrix that is needed. It's a lot of work to do when you do it by hand (don't have any maths program installed...). Will come back with my questions about that. Perhaps this should have been posted in Calculus then since many approximation methods apparently use derivatives in the process...

Please tell me if I'm on the right track so far...

 

[SevenSeas]

And you already told us:
"Now, when an idealised thread of length L isn't stretched, it will be hanging
in an arc between two points of heights H (higher height) and K (lower height)."

So does all that mean that H - K = D ?

Is there a diagram of this ETlike thing anywhere

No, D is the distance along the x axis between the poles, whereas H and K are the heights (i.e. y value) of the poles. There is a nice diagram and a lot of more explanation at the link I gave in the first post:

http://www.had2know.com/academics/catenary-equation-shape-hanging-chain.html

 

[SevenSeas]

NEAT! Thanks.

So, let u = bc and v = bD, then the equations CAN be shown this way:

e^u + e^(-u) = (H - a) / (2b)

e^(v-u) + e^(u-v) = (K - a) / (2b)

e^(v-u) - e^(u-v) + e^u - e^(-u) = L / (2b)

Agree?

Um, I think you should be multiplying, not dividing, with 2b on those right hand sides, since you switched them sides... otherwise, I think it's correct. But of course I don't understand what you're trying to do, since it tells on the page I linked that there is no algebraic solution? And most numerical methods seems to want the equations in the form

f(a, b, c) = 0 .

So then it would be:

e^u + e^(-u) - 2b(H - a) = 0

e^(v-u) + e^(u-v) - 2b(K - a) = 0

e^(v-u) - e^(u-v) + e^u - e^(-u) - 2bL = 0


I'm going to bed now, will check out thread tomorrow.

As I said, I started with trying to implement Newton-Raphsons generalised numerical method for this. It needs a matrix with the partial derivatives of these three left hand expressions, for all the three variables [i.e. a 3x3 matrix]. So I started to calculate those derivatives, using chain rule and quotient rule and stuff. Maybe there are programs that can calculate those derivatives for me, though. And then you need to find the inverse of that matrix somehow. And using that somehow, you can iterate (by hand or by computer), getting better and better results for the vector (a, b, c). But I need to know whether Newton's method is a good choice for this.

 

[SevenSeas]

I think I have an implementation of Newton's method for 3 variable functions now, but the problem is that I run up against the limited floating point precision of Javascript, so I get values like Infinity and -Infinity and then when I use those values further I get NaN (Not a Number).

The Wikipedia article about Newton's method (http://en.wikipedia.org/wiki/Newton's_method#Nonlinear_systems_of_equations) says this:

Nonlinear systems of equations

k variables, k functions

One may also use Newton's method to solve systems of k (non-linear) equations, which amounts to finding the zeroes of continuously differentiable functions F : R^k → R^k. In the formulation given above, one then has to left multiply with the inverse of the k-by-k Jacobian matrix J_F(x_n) instead of dividing by f '(x_n).
Rather than actually computing the inverse of this matrix, one can save time by solving the system of linear equations

for the unknown x_n+1 − x_n.

The last thing there is perhaps relevant here, since I think it is when I try to compute the inverse of the Jacobian matrix that numbers get too big in the process. But I don't know if a linear equation solver would do better in this respect, and whether it's any easy to implement or import.

I have tested my matrix operations separately, and they seem to be functioning correctly, although the inverse calculation suffers from limited precision so when I then multiply a matrix with its inverse there are some very small numbers in some places where there should be zero.

 

[SevenSeas]

Now I have implemented ("borrowed") an equation solver for three variable three equation systems. But when I used it, it found no solutions to my equation system. Then I discovered some errors in the derivatives in the Jacobian matrix function that I created, and also some error in the F function! But now I think I have those errors sorted out, but I still get no solutions with the equation solver!

So please check these out if they're correct (I checked the derivatives in the Jacobian matrix with an online derivative solver, but you never know).

I take it that our three equations can be rewritten on the form:

1. e^(bc) + e^(-bc) - 2Hb + 2ab = 0
2. e^(bD-bc) + e^(bc-bD) - 2Kb + 2ab = 0
3. e^(bD-bc) - e^(bc-bD) + e^(bc) - e^(-bc) - 2Lb = 0

Let F(a,b,c) be defined as

F(a,b,c) = (F1(a,b,c), F2(a,b,c), F3(a,b,c))

where
F1 = e^(bc) + e^(-bc) - 2Hb + 2ab
F2 = e^(bD-bc) + e^(bc-bD) - 2Kb + 2ab
F3 = e^(bD-bc) - e^(bc-bD) + e^(bc) - e^(-bc) - 2Lb

We are searching for a solution to the equation F(a, b, c) = (0, 0, 0)

In the process we need the "Jacobian matrix" of F, which is defined as

dF1/da	dF1/db	dF1/dc
dF2/da	dF2/db	dF2/dc
dF3/da	dF3/db	dF3/dc

As I calculate it, I get the values in J to be:

dF1/da = 2b
dF1/db = ce^(cb) -ce^(-cb) - 2H + 2a
dF1/dc = be^(bc) - be^(-bc)

dF2/da = 2b
dF2/db = (D-c)e^((D-c)b) + (c-D)e^((c-D)b) - 2K + 2a
dF2/dc = be^(bc-bD) - be^(bD-bc)

dF3/da = 0
dF3/db = (D-c)e^((D-c)b) - (c-D)e^((c-D)b) + ce^(cb) + ce^(-cb) - 2L
dF3/dc = -be^(bD-bc) - be^(bc-bD) + be^(bc) + be^(-bc)

But apparently something is wrong...

 

[Deleted member 4993]

Newton's method (solver basically uses this method) is unstable away from the solution or at J ≈ 0.

So you need to make "good" initial guess to be able to converge to the solution.

 

[SevenSeas]

Newton's method (solver basically uses this method) is unstable away from the solution or at J ≈ 0.

So you need to make "good" initial guess to be able to converge to the solution.

OK, thanks for telling me that. I found some more bugs of course, but it still doesn't converge, so that may be it. Is all my work so far in vain then, or does someone have any good idea how I might appreciate the initial values of a, b, c? The page that I linked suggests that C = D/2 when H=K, so I set that as that as init value for c. But it doesn't help. Of course there may be lots of more bugs in my code too.

Thanks for your interest.

 

[SevenSeas]

Oh my God, I have solved it!!!

It just needed more iterations to converge!

I calculated the initial values for a, b, c like this, maybe this was overkill, I don't know, haven't tested with random values yet:

c0 = D/2

b0 was calculated with Newton's method for one variable, using equation 3 (which does not have a in it but is just in terms of b and c) and putting 1 and c0 as initial values for b and c in it.

a0 was calculated as the mean of a1 and a2, where a1 was the value of a in terms of b and c (in the form of b0 and c0) in terms of equation 1, and a2 in the same way as the value of a in equation 2 in terms of b and c in the form of b0 and c0.

But now I will just have to try to speed up the calculations enormously since so many iterations were needed... I know some optimizations I can do but then we'll see.

If anyone is interested I could upload an animated gif somewhere, of the rope hanging from two moving points... but not tonight, it is late...

 

[SevenSeas]

Agree.
Let W = H + L - K
Then (screwing around with those 3 equations):
[e^(bc) - e^(b(c-D))] / b = W

So down to 2 unknowns (b and c) to be searched for.
Unknown "a" then calculated by default....using equation 1. or 2.

That's all you get from me...got nothing else

That's a really good one! As you see above, I have now solved it using approximation for 3 unknowns simultaneously. But since I need to calculate this each frame, I need something fast, and I guess your equation would be much faster to search for the answer to. So I will try to implement that instead, reusing much of the Newton's method code that I wrote.

Thanks!

But wait, do I only have one equation now, or can I use equation 3 together with your equation for b and c? (Equation 3 too is only in terms of b and c). I guess so. I need two equations for solving a system of two variables. But I guess you got that from combining (1) and (2) only, so it should be OK then. Great!

 

[SevenSeas]

For the record, here's how I ended with that:

Let u = bc, v = b(D - c) 

2b(H - a) = e^u + e^(-u) {1}
2b(K - a) = e^v + e^(-v) {2}
2b(L) = e^u - e^(-u) + e^v - e^(-v) {3}

Equating {1} and {2} in terms of 2ab:
H - K = [e^u + e^(-u) - e^v - e^(-v)] / (2b) {4}
Rearranging {3}:
L     = [e^u - e^(-u) + e^v - e^(-v)] / (2b) {3}

Adding {3} + {4} and let H-K+L = W:
Wb = e^u - e^(-v)
Substituting back in:
Wb = e^(bc) - e^[-b(D - c)]

W = {e^(bc) - e^[b(c - D)]} / b

I tested it with the help of the calculator at the site you supplied.
"Seemed" to work ok.
Hope it helps out...

Yeah, I figured that you did something like that. I actually think the step of adding {3} and {4} was unnecessary, I think you lose some information there. But your idea of eliminating a before doing approximation was genius!

Actually, today I re-implemented the Newton method for just two variables (i.e. b and c), using {3} and {4} as two separate equations in b and c, on the form:

3. e^(bD-bc) - e^(bc-bD) + e^(bc) - e^(-bc) - 2Lb = 0
4. e^(bc) + e^(-bc) - e^(bD-bc) - e^(bc-bD) + 2b(K - H) = 0

This way, the Jacobian matrix J got smaller and the linear equation solving of the

J x (x_n - x_n-1) = -F(n-1)

got much easier and faster. So I've already implemented this, and the resulting code is much faster, fast enough for my needs I think. So I consider this solved for my part.

 

[SevenSeas]

WHEN???!!!!:cool:

I'd sure like to see that in action...

Oh, sorry! The gif didn't turn out as smooth I had hoped, due to my old computer being too slow to record its own screen without hacking. But here it is (I think it will give some rough idea of it):

http://forums.tigsource.com/index.php?topic=42623.msg1052928#msg1052928