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:Say we know the state of the system (containing both grey and black boundaries) at some time $t$.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$.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 exclusionMaybe 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?
Resources exploring the problem of volume exclusion?
numerical analysis;numerics;reference request;time integration
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.