Theoretical upper force limit for ODE particle system?

Started by
3 comments, last by RandomBystander 14 years, 10 months ago
Hello everyone, I'm wondering if there is a way to actually calculate the highest possible force a particle system can hande without blowing up? I realise that this at least depends on the used integration scheme (or at least its order of magnitude), the timestep, and particle mass. Other than that I'm lost. From experience I know roughly what spring stiffnes and damping coefficients I can use and still be safe, but I never feel 100% guranteed that the simulation will not blow up. If it's actually possible to precalculate the highest allowable force, then I can clamp the upper force limit to that number and live happily ever after :-) Cheers, Mike
Advertisement
Hi Mike

Limiting the forces on your particles is probably not enough. Even at low forces things might get out of hand over time. Consider an object on ground, with slowely increasing angular velocity due to a low constant torque. The object spins faster and faster and will eventually blow up (a collision with another object will most likely cause it).

Solution 1:
if ( velocity.lenght() > THRESHOLD )
velocity *= (THRESHOLD / velocity.lenght())

Solution 2:
Introduce air friction to your particles and limiting the forces that are applied. This also gives a max velocity (terminal velocity)

Cant really answer your question, but I hope my thoughts about this helps a bit.


Hi glSmurf()

After having thought about it, I can see that you're right. It's not simply the numerical value of the force that causes instability but rather the change of force (magnitude or direction or both) over time. I'm sure there's a technical term for this, df/dt.

I'm almost certain I saw a small article on this once, but I'm not 100% sure. It was something about numerical instability of different types of integration schemes, and the author - I'm almost sure - used an equation to describe some maximum threshold or limit below which each integration scheme would behave nicely.

Anyway, as you suggest I've often used a velocity squared air friction to calm things down. It ha almost no effect at small velocities but then suddently has a huge damping effect once velocities get out of hand. The thing is, I can't always use this method, since I'm checking the quality of the physics engine via energy conservation.
Assuming it's relatively painless to calculate the force at any given time t for your particles (ie: no spring meshes), and the forces are all continuous (bouncing a particle off a wall has to be handled separately), you can use an adaptive timestep to handle any forces you like, no matter how strong.

Basically integrate with timestep t, then integrate twice with timestep t/2. If the difference is larger than some threshold, you halve the timestep and try again. If they're within some threshold, you double the timestep and try again. You should be able to arrive at a very efficient, and numerically robust, way to integrate the motion. And of course you don't have to halve/double the timestep of other particles, so the integration is the proper step size for your particles on a case by base basis.
[size=2]Darwinbots - [size=2]Artificial life simulation
Quote:Original post by h4tt3n
I'm wondering if there is a way to actually calculate the highest possible force a particle system can hande without blowing up? I realise that this at least depends on the used integration scheme (or at least its order of magnitude), the timestep, and particle mass. Other than that I'm lost.
Cheers,
Mike


My memory is a bit hazy on that matter, so all is IIRC.

In a sense, yes. Keywords for this question are for example A-stability and Dahlquist.
To answer this question for a particular integration scheme amounts to constructing a characteristic function and checking where in |C its magnitude is smaller than 1. Consider the test equation y'=-lambda.y to understand why: The exact solution converges to zero. A numerical solution must do the same: write the scheme as y(k+1)=f(y(k),dt), then this means that |f(y(k),dt)| < |y(k)|. Most numerical methods are linear in y, so one can write f(y,t)=y.R(lambda.t), where p is the function we are interested in. So, with z:=lambda.t our criterion is |R(z)|<1.

Notes: All explicit schemes (read: the simple ones) get unstable if lambda gets large enough. All explicit Runge-Kutta-methods of order s have the same R, and it's a polynomial of degree s, namely the Taylor-polynomial of exp(z). I can't find any links to actual graphs at the moment, though.

This topic is closed to new replies.

Advertisement