Nested for loop not functioning correctly for 2D particle collisions

Hello, I have very basic level knowledge of matlab and am trying to code 2D particle collisions with boundary walls and other particles. I am pretty sure that I have done the wall collisions correctly, because with just the wall collisions the graphs make sense. But when I add my particle collision code, the data gets really weird (ex: one particle has an x position of -0.0113 x 10^4 when it shouldn't be able to travel out of the boundary.)

I think it has something to do with the nested for loop, as when I adjust it, it starts giving me data that makes more sense, but I am unsure of what exactly is going wrong. If someone could take a look at this code and give suggestions as to what I could fix, that would be greatly appreciated!

(I apolgize if this is long for an ask, I just have no other way of getting feedback right now)

Initial Data (shouldn't have anything wrong here):

% number of particles
N = 3;
 
% set limits
ground = 0;
ceiling = 50;
left_wall = 0;
right_wall = 50;

collisions = 0;
c = zeros(1,1001);

% data for euler's method:
t = linspace(0,100,1001); % time from t=0 to t=100 with stepsize of 0.1s


% initial velocities for all particles
Vx = unifrnd(-10,10,N,1)

Vx = 3×1

3.1286 3.2664 -2.9203

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Vy = unifrnd(-10,10,N,1)

Vy = 3×1

-7.9536 9.1400 -6.6170

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% initial positions for all particles
xi = unifrnd(5,35,N,1);
yi = unifrnd(5,35,N,1);

x_f = zeros(N,1001);
y_f = zeros(N,1001); 
x_f(:,1) = xi;
y_f(:,1) = yi;

% c = (x1 - x2).^2 + (y1 - y2).^2; % condition for collision of p1 and p2
r1 = 8; % radius of p1
r2 = 8; % radius of p2
rad = (r1 + r2).^2; % radius of collision between any particle

Collisions:

figure
hold on;

for j = 1:numel(t)-1 % time
    for i = 1:1:N % Nth particle
    


       % wall collision for any particle
       if x_f(i,j) <= left_wall || x_f(i,j) >= right_wall
           Vx(i,1) = Vx(i,1) .* -1;
       end
       if y_f(i,j) >= ceiling || y_f(i,j) <= ground
           Vy(i,1) = Vy(i,1).*-1; 
       end
       
       % Particle-Particle Collision (PROBLEM AREA)
       for K = 1:N
           for P = K+1:N % I think this is what is going wrong!!!!!

               c(j) = (x_f(K,j) - x_f(P,j)) .^2 + (y_f(K,j) - y_f(P,j)).^2; % this part also seems a bit incorrect
               if c(j) <= rad

                   x1_coord = [x_f(K,j),y_f(K,j)];
                   x2_coord = [x_f(P,j), y_f(P,j)];

                   Vel_p1 = [Vx(K,1),Vy(K,1)];
                   Vel_p2 = [Vx(P,1), Vy(P,1)];
                   
                   % changing velocity
                   Vel_p1 = Vel_p1 - ( ( (dot( (Vel_p1 - Vel_p2),(x1_coord - x2_coord) )) / (norm(x1_coord - x2_coord)) )* (x1_coord - x2_coord) );   
                   Vel_p2 = Vel_p2 - ( ( (dot( (Vel_p2 - Vel_p2),(x2_coord - x1_coord) )) / (norm(x2_coord - x1_coord)) )* (x2_coord - x1_coord) );
                   % extracting Vx and Vy
                   Vx(K,1) = Vel_p1(1,1); Vy(K,1) = Vel_p1(1,2);
                   Vx(P,1) = Vel_p2(1,1); Vy(P,1) = Vel_p2(1,2);
                   collisions = collisions + 1;
               end
           end
       end
                   
     
    
    x_f(i,j+1) = Vx(i,1).*0.1 + xi(i,1);
    xi(i,1) = x_f(i,j+1);

    y_f(i,j+1) = Vy(i,1).*0.1 + yi(i,1);
    yi(i,1) = y_f(i,j+1); 
   end
end

What I'm using to plot:

figure(1)
for r = 1:1:N
 plot(x_f(r,:),y_f(r,:),'r')
 hold on
end

[Torsten]

Maybe you could explain in your own words what you are trying to do in the code and how you try to achieve this. It's always difficult or - if the code is errornous - even impossible to reconstruct this from the code.

[Anna]

Yes of course, thank you for answering! I have randomnly set the initial positions and velocities of an N number of particles in a set boundary, and have set two for loops which define what particle is being "observed" (for i = 1:N) and at what time (j = 1:numel(t)-1). I am using if statements to redirect the velocities if the position of each particle touches the boundary or if two particles get within a certain distance of eachother (particle collision). Once the velocity is redirected, the position of the particle changes using Euler's Method, which is at the very bottom of the overall for loop. Oh also, the positions of all of the particles are put into two different N x 1001 matrices so that I can extract that data when I need it. It's a similar thing with the velocities.

So far the only part of the code that is the problem is this part:

% Particle-Particle Collision (PROBLEM AREA)
       for K = 1:N
           for P = K+1:N % I think this is what is going wrong!!!!!
               c(j) = (x_f(K,j) - x_f(P,j)) .^2 + (y_f(K,j) - y_f(P,j)).^2; % this part also seems a bit incorrect
               if c(j) <= rad
                   x1_coord = [x_f(K,j),y_f(K,j)];
                   x2_coord = [x_f(P,j), y_f(P,j)];
                   Vel_p1 = [Vx(K,1),Vy(K,1)];
                   Vel_p2 = [Vx(P,1), Vy(P,1)];
                   
                   % changing velocity
                   Vel_p1 = Vel_p1 - ( ( (dot( (Vel_p1 - Vel_p2),(x1_coord - x2_coord) )) / (norm(x1_coord - x2_coord)) )* (x1_coord - x2_coord) );   
                   Vel_p2 = Vel_p2 - ( ( (dot( (Vel_p2 - Vel_p2),(x2_coord - x1_coord) )) / (norm(x2_coord - x1_coord)) )* (x2_coord - x1_coord) );
                   % extracting Vx and Vy
                   Vx(K,1) = Vel_p1(1,1); Vy(K,1) = Vel_p1(1,2);
                   Vx(P,1) = Vel_p2(1,1); Vy(P,1) = Vel_p2(1,2);
                   collisions = collisions + 1;

A short breakdown of this part is that the nested for loop is used to be able to compare the positions of two diff particles against each other at the same time, the c(j) equation gives the distance btw the two particles, and rad is the smallest distance particles can be apart before they collide. If they "collide" then temporary position and velocity vectors are created using the positions and velocities at the time of the collision.

Whenever I add this part of the code, the position data doesn't make sense, so obviously something is wrong. I am wondering if I don't actually need to created two new for loops, and could instead create the nested for loop at the very top of the code, or if the problem could be with the c(j), because I was never very confident in that part.

Again, thank you for answering!

[Torsten]

What is your aim in each time step ? Determine x- and y-velocities for each particle that avoid a collision with any one of the other particles ?

What if a particle will collide with more than only one particle in the next time step ? Will the adjustments

% changing velocity
Vel_p1 = Vel_p1 - ( ( (dot( (Vel_p1 - Vel_p2),(x1_coord - x2_coord) )) / (norm(x1_coord - x2_coord)) )* (x1_coord - x2_coord) );   
Vel_p2 = Vel_p2 - ( ( (dot( (Vel_p2 - Vel_p2),(x2_coord - x1_coord) )) / (norm(x2_coord - x1_coord)) )* (x2_coord - x1_coord) );                   

hinder collision with all particles satisfying c(j) <= rad ?

[Anna]

Each step is supposed to replicate a collisions, so rather than determining what x and y velocity would avoid a collision, I'd describe it as "how would the velocity change when a particle colldies with another particle?"

I should have mentioned this earlier, but I was provided a manual to help guide me. At first I had to figure out how to simulate the collisions between to particles/the walls, which I achieved. So everything I have written out was orignially written for two particles, now all I am doing is trying to adjust that code to include multiple particles.

In terms of collisions between multiple particles, I didn't actually keep that in mind while coding. Now that I am thinking of it, I think it is possible for that to have an affect on the collisions, but I am pretty sure that is not the main problem, given that from the start the data doesn't make sense. Sometimes I get NaN or 0 entries for a whole row of x/y positions for a single particle, meaning that my for loop can only detect collisions between some particles, and after a certain point it starts to malfunciton. That is what I am unsure of how to fix.

[Torsten]

At least the loop that detects wall collisions for the particles should be placed after the double loop that detects particle-particle collisions. Otherwise, necessary reversions in direction are simply overwritten in the particle-particle collision loop.