Rope Physics Help

[NotARaptor]

ZIP files are case sensitive. You told it to require "Rope". But the file is named "rope". Heads up to anyone downloading this to unzip it and change that code to make it run.

Sorry about that! Very new to LUA, should have thought about that. I keep indexing arrays from zero too!

[zorg]

Sorry about that! Very new to LUA, should have thought about that. I keep indexing arrays from zero too!

(And it's also called lua, not LUA; it's not an acronym, it's portuguese for moon. )

[NotARaptor]

(And it's also called lua, not LUA; it's not an acronym, it's portuguese for moon. )

Noted Thanks.

@pgimeno - I did start composing a reply about different integrators, timesteps, precision issues, spring stiffness and why the relaxation technique used with Verlet is so popular... but it ended up sounding really preachy and boring, so I didn't post it. If you do need any help with it let me know

[pgimeno]

I'd appreciate if you post it I love to learn.

[NotARaptor]

Ok... <cracks knuckles> let's do this.

Boiling the problem down to the basics, we have a system of N objects - we know their positions and velocities, and their masses, at a specific point in time.

We want to simulate the system over dt time from that point, and see where they end up.

We use Newton's laws - specifically the second law here.

If we can work out all the forces acting on the system over the time dt, that force function will give us the acceleration. Since the acceleration is the derivative of the velocity with respect to time, we can integrate that acceleration function to give us the velocity. That can then be integrated over time to give us the position. Simple!

For the moment we're going to treat each object in the system individually, and just treat them as a point mass for now.

We have F(t)

F(t) = m * a [Newton's 2nd Law of Motion]

from that we have a(t)

Integral of a(t) over [0, dt] will give us v(t) [we have the initial velocity - that's the constant C your maths teacher kept going on about]
Integral of v(t) over [0, dt] will give us p(t) [we have the initial position]

Seems simple, right? WRONG

It's perfectly fine if we just consider a constant force - like (simplified) gravity. Let's do that - in this world gravity is a constant force in the +y direction, with a magnitude of G.

We'll call the initial positions and velocities p0 and v0. The t variable will be considered to vary between 0 and dt.

This means that

Force calculation
F(t) = <0, G*m> [Just using 2d here - 3d is exactly the same]

Newton's 2nd law to get the acceleration - mass cancels out
a(t) = <0, G>

Integrate to get velocity [not forgetting the constant!!]
v(t) = v0 + <0, t * G>

Integrate again to get the position
p(t) = p0 + <0, 0.5 * t * t * G>

Oh that's nice, we can find the end point!!

New velocity:

v' = v0 + <0, dt * G/m>

New position:

p' = p0 + <0, 0.5 * dt * dt * G/m>

How simple! How easy! It's perfect. Except it's kinda pointless.

See we said here : "if we just consider a constant force"

Well, most forces aren't. Most forces will depend on the object's position ([inverse squared] distance from a point, for example) or the object's velocity (drag, friction).

If the force depends on the position, the acceleration depends on the force, the velocity depends on the acceleration, and then the position depends on the velocity, then instead of a neat 1-2-3, we have a system of differential equations.

One way to sort it out is to approximate the system - if dt is "very small", then over that time we can consider that the force will be constant. If the force is constant, then the acceleration is constant, and we're back to our neat 1-2-3.

This is forwards Euler integration. We want to evolve the system over dt, so we calculate the force at t=0, that gives us the constant acceleration, and from there we compute the velocity and position like we did above.

Of course since the force won't actually be constant across dt, this will produce an error. The smaller the timestep, the smaller the error.

Image from wikipedia's article on Euler integration showing the error:



Now, how small do you need the timestep to be to minimise error? It depends - a lot. But just doing dt in one step probably won't cut it for all but the simplest simulations, you'll have to do dt in increments - for one timestep in the world (love.update call), you might have to split dt into 10, 100, 1000 or more intervals and simulate them.

And then of course you run into numerical precision issues, which I'll get into a bit more later. Basically floating point numbers are very accurate for a certain range, but very small (close to zero) and very large (far away from zero) numbers aren't represented as precisely. And when you multiply very teeny tiny numbers with very huge numbers, the result will be wrong. So you use a smaller timestep to get more accuracy in the physics, but the tiny timestep causes precision issues so you get less accuracy. It's fun.

You can also minimise the error by using a different integrator. I'm not going to go into huge detail on this here (google "RK4"), but a simple one to explain is the midpoint method. We do Euler like above, but instead of calculating the position at dt, we calculate it at half-way (dt/2), then recompute the forces and acceleration at THAT point, backtrack to the original point and use the revised computation over the whole dt interval. The error with Midpoint goes down much faster than with Euler, and even faster with higher-order methods like RK4. But...

Springs, constraints and noninterpenetration

It's fairly easy to add springs into the system we've described above. If there is a spring between two point-masses A and B, then at each point in time, the spring exerts a force both on A and B that either pulls them together or pushes them away, based on the distances between them, the rest distance of the spring, and the spring coefficient (basic Hooke's law). This is perfect example of how a force depends on the position.

This could be used to simulate a rope. The issue is with the spring coefficient - if you keep it low, it will behave like a spring - bouncy. Increase the spring coefficient to make a stiffer spring, make it high enough and it will behave like a rigid joint. Well, in theory anyway. In reality it's a bit different.

When the spring coefficient is really high, the forces generated are very strong. And vary a lot over a small timestep, when you have multiple springs push-pulling particles all over the place. This makes the numerical integration harder to calculate. MUCH harder. So... reduce the timestep? Can do, but then very small numbers multiplied by very big numbers... BOOM numerical precision issues. Like I said, fun fun fun.

Let's consider noninterpenetration. We were dealing with point masses, but if we add a radius to their description, other forces (gravity etc) won't change for them and they'll behave the same way. So now we have little discs (or balls, in 3d).

Noninterpenetration is basically "don't let them overlap". Circle-circle collision detection is easy to calculate (distance between the centres is less than the sum of the radii => intersection) and likewise collision response is easy (push them along the vector formed by the centres). But how do we do this in simulation?

In physics classes, you draw the diagram of the two balls just touching at the moment of impact, and then use the laws of conservation of momentum and conservation of energy to work out the impulses (immediate changes in velocity) that happen at that point. Then you have the velocities before the collision and the velocities after the collision.

Over the timestep dt, two discs are moving, and might intersect during that time period. Let's say we assume they do - to accurately resolve the collision, we need to determine the PRECISE time that they collide, move the objects to that position (and velocity) resolve the collision as above, then carry on simulating for the rest of the interval - where another collision might well occur.

How do we determine the exact time they collide? Well, it's easy to check if they are colliding or not at any point, so we could step through dt to find a time t1 before they collide, a time t2 after they collide, then do a recursive subdivision (like binary search) to determine the actual collision time.

We could also solve it - we have the position of both objects over the interval represented as a quadratic function of time. Given that, we can work out their distance as a function of time, and solve for distance = zero.
If we perform a change of reference so that the first object is stationary, we reduce it down even further. And we don't need the distance (nasty square root), we just need the squared distance - the distance will be zero if and only if the squared distance is zero. So that's solving this equation:

(a*t^2+b*t+c)^2+(p*t^2+q*t+r)^2=0 for t

A quick trip to wolfram alpha gives me this...

equation.png (54.42 KiB) Viewed 5469 times

My intuition says this isn't going to be a fast calculation to run multiple times per frame.

Of course because we're not dealing with true real numbers in the mathematical sense, we'll only ever have an approximation of the time. We'll have one time t1 where the discs are NOT colliding, and then another time t2 where they are JUST colliding (overlapping). Given those times, we can calculate the collision response above (how far do we have to move the discs so they don't collide), and then based on that amount we can perform a best-guess interpolation between t1 and t2 to get our (I changed this idea to the binary search a few paragraphs above - this would work too though)

With two objects this is fairly simple, but with N objects it gets harder. The common method used in physics engines is to partition the objects into "islands" of potential collisions. Maybe out of 100 objects you'll have 95 which are far away from everything else and which won't collide - these can be resolved using the normal integration step. The other 5 might be in two "islands" - one island containing A,B and C which might potentially collide during dt, and D and E which might collide with each other. You then resolve those islands as independant subsystems using the step system above. Don't forget it's totally possible for more than one collision to happen during the timestep, and they all have to be resolved EDIT: resolved in order, as well!. To be TRULY accurate you need to handle the possibility that more than two objects might collide at the same time, which makes the resolution step even higher, but in practice you can just handle collisions between two objects.

But can't we just represent the noninterpenetration constraint as really stiff springs that only act when discs are really close to each other? Of course we can! But nice stiff springs give unstable equations... so we need smaller timesteps (more calculation, slower)... but that gives us tiny numbers multiplied by huge numbers again, which will cause numerical instability. Bloody floating point!!!

Of course that's strictly for discs - to handle objects with different geometries, you have to take into account torque, rotational velocity and angle as well as linear position, velocity and acceleration. And the actual collision detection is a lot harder. AND the collision response is even harder, if you need that. Oh and that's just for 2d - going to 3d with balls instead of discs is easy-peasy, but more complicated geometries are much harder in 3d.

------------------- INTERLUDE -------------------

I've veered a bit off "rope simulation" here and gone into general simplified Newtonian physics simulation - sorry about that - but it gives a good background. The main point I'm making is that if you use classical numerical integration for simulation, in order to get any accuracy at all with systems beyond "two balls and one spring with moderate stiffness", you'll be doing a lot of calculations every timestep. And you can't just write a "generic" system that will handle anything you throw at it - you'll have to fine-tune error levels, timesteps and everything else for your problem space. If you solve it for the general case you're basically just writing your own love.physics engine.

----------------- END INTERLUDE -----------------

Enter Verlet integration with constraint relaxation!!!

Verlet integration isn't anything really special - the main point about it is that the velocity isn't calculated (or stored) directly; it's implicitly represented by the difference of the positions. This makes some things easier, some things harder.
But the reason it's used so much in game simulation is because it works well with the constraint relaxation process. The forces are all calculated in the same way as above, but noninterpenetration and other constraints are handled in a much more programmer-like fashion - we just enforce them directly on the positions. This can (read: will) cause other constraints to be violated, so we just repeat the process a few times - and like magic, it converges to a solution! Really fast too.

This can even do things like IK with practically zero effort. Let's assume we've got a simple arm set up in a similar way to the rope from the .love earlier. However many segments as you like. If you enforce as a constraint that the last point must be at a certain fixed position, on the first iteration it will do that, which will make the last segment (for example) really long. In the next iteration, that segment will enforce its length constraint on its neighbours, moving the previous point and the end point. The end point will then be moved back to the target position. As this iterates again and again, the error is "spread out" until everyone is happy (if that's possible).

By contrast, try to imagine how you'd handle that with classical numerical integration. The target point attracts the end effector with a force? You'd have no way of knowing how long it would take to get there.

The Verlet integration with constraint relaxation technique is not physically accurate - the system will lose energy over time, and momentum might not be conserved. Your physics teacher would probably disapprove (I think mine rolled over in his grave when I didn't bother with a mass variable and just multiplied Newtons by time in my other code), but it looks "good enough" for games.




To sum up, I don't think changing your integrator from Euler (or Verlet) to RK4 will help you for a rope simulation, not if you want it to run in realtime. Whichever way you choose to go, you're going to have lots of calculations to do, and the "iterating 20 times" method of Verlet+Relax is by far the least inexpensive method, even if it looks on the surface like it's inefficient. It also adds tons of benefits for collision and other constraints which are harder in the more "classical" way of doing things.



As a footnote, I should add that there are integrators that work well with stiff equations - implicit methods. Backwards Euler for example. As far as I know these aren't used in "casual" programming, because you'd need something like Matlab to solve a giant sparse matrix of coefficients for a big system of equations. But it's not something I've really looked into other than the absolute basics.



... and breathe....

Hope this is helpful! Any questions (or mistakes I've made), just let me know.

Edit #2 - fixed a rookie error with the mass calculation

[NotARaptor]

I'll make some code examples of the various steps at some point!

[Ref]

I always liked Micha's rope simulation.
Contains a lot of interesting code.
Best

[NotARaptor]

Nice one @Ref - that's another example of Verlet integration with constraint relaxation, but obviously written to specifically handle the rope and nothing else. It looks nice

Personally I'd tweak the variables a bit as it feels "floaty", but it doesn't feel "stretchy" which it would if it was written with springs.

The key-based movement of the end of the rope is done by simply setting the endpoints (one of the advantages of this method), but as you tweak the number of iterations, the "power" of the push varies, because that's set before the relaxation loop happens. Of course if the number of iterations is fixed, this doesn't matter at all.

[pgimeno]

Thanks a lot for the detailed explanation. I have a better understanding of the issues now.

If the force depends on the position, the acceleration depends on the force, the velocity depends on the acceleration, and then the position depends on the velocity, then instead of a neat 1-2-3, we have a system of differential equations.

This made many things click into place, and made me realize why my OP's question:

Can someone help me with the force calculation?

was a fool's errand. Incidentally, by pure chance I also found this Wikipedia page a few days ago: https://en.wikipedia.org/wiki/Wave_equa ... _dimension (while reading this thread) and noticed an alarming similarity between the 1D setup explained there and the 2D rope I was trying to make, which made me suspect that if so many clever people were throwing so much effort into this and having so much trouble with it, then the answer to my problem was probably not that simple. That paragraph has helped me a lot in understanding why. And also why I was getting no results no matter what I tried.

I don't need to handle collisions for my use case. Still, your help on that part may prove quite helpful in future.

Whichever way you choose to go, you're going to have lots of calculations to do, and the "iterating 20 times" method of Verlet+Relax is by far the least inexpensive method, even if it looks on the surface like it's inefficient.

You've convinced me of that, thanks. My concern is that I'm targeting lower-end computers and the rope calculation is just one part of the frame time; the game is not just about ropes, there's a lot more to process. If the rope takes, say, 1/30 of the frame time, I consider that a lot.

Next I'll experiment with finding a compromise between spring stiffness, iteration count and distance correction, and see where it leads me. So far, I've noticed that performing distance correction after 20 iterations, the result was still fairly credible. With 1 iteration, it was far too incredible because it wiggled too much and too fast.

Hope this is helpful! Any questions (or mistakes I've made), just let me know.

It was! And allow me to take your word on that...

And then of course you run into numerical precision issues, which I'll get into a bit more later. Basically floating point numbers are very accurate for a certain range, but very small (close to zero) and very large (far away from zero) numbers aren't represented as precisely. And when you multiply very teeny tiny numbers with very huge numbers, the result will be wrong. So you use a smaller timestep to get more accuracy in the physics, but the tiny timestep causes precision issues so you get less accuracy. It's fun.

I'm quite acquainted with the precision issues of floating point calculations, and I have to disagree with the part I've highlighted in boldface. No matter the magnitude of each number, to multiply two floating-point numbers, you multiply the mantissas and add the exponents, and assuming no overflow or underflow occurs, you get a decent result (granted, multiplying two n-bit numbers gives a 2n-bit number, of which you're rejecting half, but believe me, that's not that big a loss typically).

No, the problem comes with adding big and small numbers. The numerical precision problems you mention still hold, because integration is basically a sum of many tiny parts into a big part, and when the big part is big enough, the tiny parts may become too insignificant to even have an effect on the sum. For example, if your mantissa has 8 digits and you attempt to add 2 + 0.0000000000002, the exact result is 2.0000000000002 but due to the precision of the mantissa, you just get 2.0000000 which is still 2.

EDIT: The part I've marked in blue is also incorrect; floating point numbers are notable for that. You need to go to REAAAAALY low numbers (under 10^300, to the precision limit) to start to see precision loss.

I'm well aware of these problems. That was one of several reasons for my choice of representing the rope as a chain of vectors that make up the segments, as these are smaller and hopefully a wee bit closer to the magnitude of the numbers we're using (if I have 16 segments, the net gain could be up to 4 bits of extra precision). But now that I better understand the issues, I notice that such a representation is not very appropriate for calculating the spring forces in both directions.

As a final note, I took a look at this paper: http://www.kuffner.org/james/software/d ... Thesis.pdf which explains the physics of e.g. robot arms with prismatic or revolute joints, making a construction that resembles a lot the one I made with Box2D earlier in the thread. But a first quick look revealed that it probably falls into the "you'd need something like Matlab to solve a giant sparse matrix of coefficients for a big system of equations" category that you mention.



@Ref, Micha's thread is linked earlier. It uses 40 passes over all segments of the rope per frame, which I was trying to avoid.

[NotARaptor]

I'm quite acquainted with the precision issues of floating point calculations, and I have to disagree with the part I've highlighted in boldface. No matter the magnitude of each number, to multiply two floating-point numbers, you multiply the mantissas and add the exponents, and assuming no overflow or underflow occurs, you get a decent result (granted, multiplying two n-bit numbers gives a 2n-bit number, of which you're rejecting half, but believe me, that's not that big a loss typically).

No, the problem comes with adding big and small numbers. The numerical precision problems you mention still hold, because integration is basically a sum of many tiny parts into a big part, and when the big part is big enough, the tiny parts may become too insignificant to even have an effect on the sum. For example, if your mantissa has 8 digits and you attempt to add 2 + 0.0000000000002, the exact result is 2.0000000000002 but due to the precision of the mantissa, you just get 2.0000000 which is still 2.

That's... very true, makes sense I was just going on what I remember from doing it - found instability, checked the numbers and the forces were very large and the timesteps were very low. Something about that combo made the resulting maths unstable! (This was over 15 years ago though... my memory isn't perfect!)

Thank you for the more correct explanation of where that error comes from

#Edit : Thanks for the link! I will give that paper a good read later