# Resources exploring the problem of “volume exclusion”?

Consider the following situation:

There are two boundaries -- one is denoted using grey lines, and the other is denoted using black lines. The boundaries are numerically represented using "vertices", with edges assumed to be connecting these vertices, but otherwise not necessary for the calculation.

At a given vertex (say, vertex $i$), we can write ODEs of the form:

\begin{align} \frac{\textrm{d}x_i}{\textrm{d}t} &= A_i(...) \\ \frac{\textrm{d}y_i}{\textrm{d}t} &= B_i(...) \end{align}

Here, the $...$ are placeholders for various variables that could be, for example:

• the position of other vertices on the same boundary (relevant for instance, if we are modelling "spring-like" connections between vertices)
• the position of vertices on the other boundary (relevant for instance, if we are modelling repulsive forces between the boundaries, as shown in the situation diagram above as an example)
• various other variables (for example, maybe a mystery chemical signal which determines the strength of the forces acting on the vertex in some particular direction)

In order to model the dynamics of the boundaries, one can integrate the vertex ODEs. For a simple example as shown in the diagram above, we can simply consider the ODEs for all the vertices together in one large system, and then integrate this system of ODEs numerically using integrator.

## A possibly invalid integration simplification?

If the boundaries rarely ever interact (i.e. they are rarely so close as to actually affect each other, depending on how the model is set up), then one could attempt another approximation, in order to lessen the computational burden on the integrator:

1. Say we know the state of the system (containing both grey and black boundaries) at some time $t$.
2. In order to compute the state of the system at $t + \Delta t$, pick one boundary at random, call it boundary $P$ for "picked". Evolve the chosen boundary from time $t$ to $t + \Delta t$ using numerical integration, but consider the other state of the other boundary (call it $Q$ for "not picked", or "not $P$") to be held fixed during this calculation. So, when we need to use information regarding $Q$ to evolve the state of $P$, we only use information about the state of $Q$ at $t$.
3. Evolve state of boundary $Q$ via integration from time $t$ to $\Delta t$, and if any information regarding $P$ is needed for the calculation, only use information regarding $P$ from time $t$.

Summarizing, this method aims to reduce the "burden" on the integrator by integrating small pieces of the global system at a time, while assuming that the rest of the system has not evolved.

## The problem of "volume" exclusion

Maybe the boundaries denote special areas that never intersect: for instance, both boundaries might represent separate membranes which can come very close to one another, but do not intersect. I call this "volume exclusion" (two "volumes" cannot intersect), but perhaps there is better terminology.

In any case, one might implement volume exclusion by assigning repulsive forces to the vertices representing each boundary, where these repulsive forces depend on the distance between a particular vertex and the closest point away from it on the other boundary (see the situation diagrammed above for an example).

There are several issues with implementing such a rule however, with one in particular being how "intersections" (due to numerical error) are handled. If one vertex does end up in the volume contained by another boundary, how do we correct for it without breaking the integration process?

Are there any resources which explore such problems?

• How is what you're describing different from the usual explicit Euler method? I think these are called non-penetration constraints. There are ODE solvers that support event location (SUNDIALS, for example). Introducing repulsive forces may make your ODE system stiff, making the explicit Euler method unsuitable. How familiar are you with these things? – Kirill Jul 4 '15 at 6:19
• @Kirill I think what I described didn't make any assumptions regarding the type of the integrator: explicit, implicit, or whatever. Thank you so much for non-penetration constraints keyword -- that leads to some pretty good reading on the internet! Maybe you can make an answer mentioning that keyword? I would be happy to accept it as an answer. I am familiar with stiffness issues, and am using an integrator that is implicit and adaptive in method (scipy.odeint): that is, it is supposed to switch to a method better suited for stiff problems when stiffness is encountered. – user89 Jul 4 '15 at 18:41
• @Kirill However, when I turn on message reporting regarding whether a method switch decision was made, I find that no messages are printed...anyway...that's too technical and specific an issue for me to make a question out of on scicomp.stackexchange. – user89 Jul 4 '15 at 18:42

1. I believe the usual name for these things is non-penetration constraints. As you say, there is already quite some literature on this, especially in robotics and computer graphics.

2. Your method looked like explicit Euler to me because I didn't notice that you didn't say how you integrated the non-fixed half of the system. So you seem to know what you're doing.

Regarding the question of whether it's valid, here is a simplified analysis of it. Take the coupled ODE $$\dot u = f(u,v), \qquad \dot v = g(u,v).$$ With a time step $\delta t$, the exact solution at $t+\delta t$ is $$u(t+\delta t) = u(t) + \int_t^{t+\delta t} f\big(u(s),v(s)\big)\,\mathrm{d}s.$$ The solution with "fixed" $v$ is $$\tilde u(t+\delta t) = u(t) + \int_t^{t+\delta t} f\big(u(s), v(t)\big)\,\mathrm{d}s,$$ which is approximated by the corresponding step of the ODE solver of order $s$: $$\tilde u(t+\delta t) = u^{n+1} + O(\delta t^{s+1}).$$

The difference between the two schemes is: $$u(t+\delta t) - u^{n+1} = O(\delta t^{s+1}) + \int_t^{t+\delta t}\Big( f(u(s),v(s)) - f(u(s),v(t))\Big)\,\mathrm{d}s,$$ and the second error term here can be approximately bounded using Taylor series: $$O(\delta t^{s+1}) + |f_v||\dot v|\delta t^2.$$

The point is that the scheme is first-order accurate, and, depending on how small $|f_v|$ is (i.e., whether the interaction between the two system parts is the dominant interaction), it might be sufficiently accurate for your purposes. It might break down for a sufficiently stiff ODE, but that's a separate question.

3. For the purposes of looking up more references, I think your method is very close to a first-order split-step method. In a typical split step method, if you have a linear system $$\dot y = \mathbf{A}y, \qquad y(t) = e^{t\mathbf{A}}y(0),$$ and you define a splitting $\mathbf{A} = \mathbf{A}_1 + \mathbf{A}_2$, you can find appropriate constants so that $$y(t) \approx e^{c_1 t\mathbf{A}_1}e^{c_2 t\mathbf{A}_2}\cdots y(0).$$

For the linearized system $$\mathbf{A} = \begin{pmatrix} \mathbf{A}_{11}&\mathbf{A}_{12}\\ \mathbf{A}_{21}&\mathbf{A}_{22} \end{pmatrix},$$ your method is close to being the split-step method with the splitting $$\mathbf{A} = \begin{pmatrix} \mathbf{A}_{11}&\mathbf{A}_{12}\\0&0\end{pmatrix} + \begin{pmatrix}0&0\\ \mathbf{A}_{21}&\mathbf{A}_{22}\end{pmatrix}.$$ Notice that both $e^{t\mathbf{A}_1}$ and $e^{t\mathbf{A}_2}$ evolve only part of the variables, and keep the rest fixed. The splitting operators are singular here, which is a bit unusual. This isn't quite the same as your method because $e^{h \mathbf{A}_2}e^{h \mathbf{A}_1}$ would use the new value of the part of the system that gets evolved first, rather than the old value.

Formulating it this way would make it easier to analyze the method's order and stability properties. For your purposes the first-order split-step method might be sufficient, in which case formulating it in the generic terms of split-step methods might be a bit overkill.