I want to numerically simulate the motion of a particle (unit mass) in $N$-dimensional space according to the force $f(\vec{x})$ and bounded by a set of linear constraints (equalities as well as inequalities):

$$A_i \vec{x} = \vec{c_i}$$

$$B_j \vec{x} \leq \vec{d_j}$$

The behavior of the particle when it comes in contact with a "wall" formed by an inequality is a perfectly inelastic collision.

What I'm thinking right now is something like this:

for step in range(total_time//dt):
    p += f(x) * dt
    x += p * dt
    if violate_constraints(x):
        x, p = correct(x, p)

It is easy to check if $\vec{x}$ violates any constraints, but I'm not sure how to implement the correct(x, p) part.

I think the correct behavior for correct(x, p) is to put the position $\vec{x}$ back to the closest point in the "allowed" space, and project $\vec{p}$ onto the intersection of subspaces according to which constraints $\vec{x}$, before correction, broke. However,

  1. I have no proof if this is correct
  2. I don't know how to put $\vec{x}$ back to the closet point in the allowed space

EDIT: Putting $\vec{x}$ back to the closet point in the allowed space is what I think should correctly handle the "corner" cases where the particle comes in contact with 2 or more walls.

  • 1
    $\begingroup$ "inelastic" do you mean elastic? In a perfectly inelastic collision with the wall the projectile sticks to the wall and stays there. $\endgroup$ Sep 13, 2022 at 12:25
  • $\begingroup$ @SeanE.Lake Yes I mean inelastic. I want the particle to stop its motion in the normal direction of the wall when they come in contact. That said, I think the difficulty of this problem comes from the fact that I want a general solution for any orientations and positions of the walls. With minor modifications, a method that works for elastic walls should work for inelastic walls as well. $\endgroup$
    – haoyu
    Sep 13, 2022 at 16:49
  • 1
    $\begingroup$ I’m voting to close this question because it belongs to the Computational Science site. $\endgroup$
    – Miyase
    Sep 13, 2022 at 19:27
  • $\begingroup$ In a perfectly inelastic collision friction kicks in, so all motion would stop on a collision with a wall. $\endgroup$ Sep 13, 2022 at 22:12
  • $\begingroup$ @SeanE.Lake Aren't collisions perpendicular interactions? Friction as a parallel interaction can be treated separately, right? $\endgroup$
    – haoyu
    Sep 14, 2022 at 6:26

2 Answers 2


Let's say you write your vectors as lists with length n. Then I would write something like this

x=[...] #your initial state
def f(x):
    f = []
    return f

def evolve_to_next_state(x, p):
    previous_x = x
    previous_p = p
    for d in range(len(x)):
        p[d] += f(x)[d]*dt
        x[d] += p[d]*dt
    for d in range(len(x)):
        if violate_constraints(x[d]):
            x[d], p[d] = correct(d, previous_x, previous_p)

 def correct(d, previous_x, previous_p):
    x[d] = previous_x[d]
    p[d] = -previous_p[d]
    return x[d], p[d]

This should solve the corner problem, going through each dimension singularly. Don't forget that now violate_constraints takes floats as input. You then have to write a function that plots it, or at least goes through each time interval... Hope it could help

Edit: 1. assume point particle... 2. Obviously this corrects only the inequality, you have to implement the code to correct also the linear equality condition

  • $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Sep 13, 2022 at 19:19

With a perfectly inelastic collision between the wall and the mass, the mass will stick to the wall at the point of the collision and all of its velocity drops to 0. That part of the problem should be straight forward to implement (e.g., if collision_detected: v = 0). The goal is then to find the point at which it collided with the wall and set its position to that point.

Using the image below, produced by Nikos M on this StackOverflow post, we can see that if $x_\text{mass}>x_\text{wall}$, then we can compute the depth, $D$, past the wall's edge via, $$ D=x_\text{mass}-x_\text{wall}+r. $$ where I assume the mass is a ball of radius $r$ (obviously if it's a point mass, then $r\to0$ is simple enough a task).

You still have your previous position increment, $p\cdot\mathrm{d}t$, so your $x$ distance is reduced from $p_x\cdot\mathrm{d}t$ to $p_x\cdot\mathrm{d}t-D$. Simple geometry should lead you to the angle of attack, $\tan\theta=p_y/p_x$, so then you can determine the height you can go, $$H=p_y\cdot\mathrm{d}t-p_x\cdot\mathrm{d}t\tan\theta$$ Thus, you have your new positions when crossing the wall along the $x$ axis: $\vec{x}_\text{new}=D\hat{\mathrm{e}}_x+H\hat{\mathrm{e}}_y$.

enter image description here

The case for crossing the $y$ axis is done with the same procedure. For the corner cases, it seems to me that you can simply place the mass at the corner itself, that is, $$\vec{x}_\text{new}=(x_\text{wall}-r)\hat{\mathrm{e}}_x+(y_\text{wall}-r)\hat{\mathrm{e}}_y,$$ as the collision with the two walls would place it directly in the corner, rather than trying to do a sweep in the $x$ and then $y$ direction (though that might be necessary for cases where the collision is near-enough to the corner, I'm not sure).

  • $\begingroup$ Nope. What you want to do is, "If violates inequalities, find the vector that points from the boundary to the object and subtract twice that vector from the position." That correctly handles the collision no matter when in the intermediate time it happened (I'm not sure if there are problems in the literal corner cases). That still leaves the equation constraints, but those should be satisfied by a coordinate transform. $\endgroup$ Sep 13, 2022 at 12:45
  • $\begingroup$ What if there are multiple walls in multiple different directions? Like @SeanE.Lake commented, I don't know how to correctly handle corner cases. Also, the walls in my problem are perfectly inelastic. $\endgroup$
    – haoyu
    Sep 13, 2022 at 17:02
  • $\begingroup$ @SeanE.Lake ah, yes. that is a good point. Will work on an edit now (was at work the last 11 hours) $\endgroup$
    – Kyle Kanos
    Sep 13, 2022 at 23:59
  • $\begingroup$ @haoyu I somehow overlooked the inelastic part of the query. I'll make an edit in that regard. $\endgroup$
    – Kyle Kanos
    Sep 14, 2022 at 0:00
  • $\begingroup$ @haoyu Why bother with a simulation, then? They're all just going to adsorb into the corner their initial velocity vector is most pointed at (if they hit had slide, that is, instead of hit and stick). $\endgroup$ Sep 14, 2022 at 0:43

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