2
$\begingroup$

I have a system of ODEs for functions $f_i(t)$. At each time $t$, $f_i(t)$ is the position of particle $i$. The functions $f_i$ have a monotonicity property: at all times $0 < f_1(t) < f_2(t) < \dots < f_n(t)$. This property is preserved by the true dynamics, since there is a divergent repulsive term in the ODEs when any two $f_i$ get close together. But they do, in fact, get close together. In fact they get close enough together that crossings can easily happen in a numerical simulation (and if they do then the rest of the run is ruined).

I have been simulating my system with ode15s in MATLAB (all the other solvers fail quite spectacularly). My method has been to use an events function to search for collisions, and if I find a collision, then I just let the simulation fail. (This is reasonably convenient to program, because I need the events function anyway for an unrelated reason.) I am not sure how efficient this approach is. Would it be more efficient to construct ODEs for $g_i=f_i-f_{i-1}$ (with the convention $f_0=0$) and then utilize the NonNegative option, then just reconstruct $f_i$? Or would that end up being implemented in basically the same way internally? (I don't really know much about how odezero works.)

Improve this question
$\endgroup$
11
  • $\begingroup$ Does it actually describe a mechanical system? Do f_i represent individual particles? Or distribution functions? $\endgroup$
    – Maxim Umansky
    Sep 11 '20 at 15:37
  • $\begingroup$ @MaximUmansky They are particles, each with a defined position at each time. I edited the first paragraph to make this more explicit. $\endgroup$
    – Ian
    Sep 11 '20 at 15:40
  • $\begingroup$ In a sense, a divergent repulsive force is the same as a collision. By adjusting the form of the force you can keep the [centers of] particles at desired distance apart. $\endgroup$
    – Maxim Umansky
    Sep 11 '20 at 15:58
  • $\begingroup$ @MaximUmansky mean that it diverges as the distance goes to zero, very similar to (but not quite the same as) the usual electrostatic repulsion between same-charged particles. This at least should cause there to never be a collision. $\endgroup$
    – Ian
    Sep 11 '20 at 16:03
  • 1
    $\begingroup$ It looks like the problem is "overshooting" which allows the trajectories to get close to each other in spite of the [stiff] repulsive force term. Probably using an implicit integrator like ode15i would prevent it. $\endgroup$
    – Maxim Umansky
    Sep 11 '20 at 18:51
2
$\begingroup$

One way to approach this is to use some adaptive scheme for your timestepping. If you have a fixed timestep width, it can happen that two of your particles switch positions (as you said). If you build a mechanism into your algorithm that reduces the time-step $h$ whenever two particles get too close to each other $|f_i(x) - f_{i+1}(x)|<\epsilon$, then your divergent force will prevent that they infact switch places. Depending on your timestepping scheme that should not be too hard to implement.

You could use some heuristic like:

$h = min(h_0~~ , ~ \omega \min{\Big( |f_i(x) - f_{i+1}|\Big)})$

and play around with the parameter $\omega$ until it runs sufficiently stable. With this you also have an upper bound of your timestepping error.

Another approach is to test for the correct ordering after each timestep. If you find that two particles have changed places, then you discard the last step and redo it with a timestep that is a lot smaller.

I should also mention that the two approaches will not scale well if the number of particles gets too high, as then the likelihood of any two particles being close to each other at each instance of time is quite high.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Maybe if you browse through the literature for Molecular dynamics (or their numerical methods), you might find how they overcome this problem in higher dimensional simulations. As far as I understand, any 2D/3D MD simulation with repulsive forces should have the same set of problems you describe. But it is also entirely possible that they simply ignore these events if they are interested in some statistical quantity. What are you trying to measure?:-) $\endgroup$
    – MPIchael
    Sep 14 '20 at 8:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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