How Do I Update Particle Positions in a Constrained 3D Simulation?

[logix88]

I am trying to simulate a 3D particle simulation, where there are a 1000 particles of random x,y,z postn varying area between 2 n 4, and mass... There are different forces acting on each particles, Force due to heat, wind, and gravity, we calc acceleration from that, after which change in distance and hence new x,y,z positions are obtained. I have written most of it (at least I thnk I did),

I need help updating these new x,y,z postion at time range (1:10), I am sort of stuck there.

I also need to lock down this 3D space to a 100x100m box, it seems to go beserk,...

and finally, these particles should get stuck to the top and bottom if they reach there... and if it goes beyond left or right, it should appear in the opposite side.

Here's the code, let me knw what you guys think.

NB : particle data file here : http://bit.ly/6qJaNb

%n = no.of particles  REPLACE 1000 with n

data=load('particledata.dat');

scatter3(data(:,1),data(:,2),data(:,3),data(:,5),'filled'), view(-60,60)

m=data(:,4);

g=9.81;

x=data(:,1);
y=data(:,2);
z=data(:,3);

F_w=zeros(3,1,1000);

F_g=zeros(3,1,1000);

F_h=zeros(3,1,1000);q=0.08;

A=data(:,5);            
            
for t=2:1:10
    for w=1:1:1000;
     
            F_g_z(w)=-m(w)*g*z(w);   % force due to grav

            F_g(1,1,w)=0;
            F_g(2,1,w)=0;
            F_g(3,1,w)=F_g_z(w);

            rx=50-y(w);
            ry=x(w)-50;

            F_w_x(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*rx;    % force wind
            F_w_y(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*ry;
            F_w(1,1,w)=F_w_x(w);
            F_w(2,1,w)=F_w_y(w);
            F_w(3,1,w)=0;

            T(w)=1/16*(70*x(w)+(50-y(w))^2);

            F_h1(w)=q*A(w)*T(w)*z(w);             %force heat

            F_h(1,1,w)=0;
            F_h(2,1,w)=0;
            F_h(3,1,w)=F_h1(w);

            F=F_g+F_w+F_h;
        
            a=F/m(w);                     % acceleration
            
            
            u(1,1,w,1)=a(1,1,w);
            u(2,1,w,1)=a(2,1,w);
            u(3,1,w,1)=a(3,1,w);
            
            u(1,1,w,t)=u(1,1,w,t-1)+a(1,1,t)*t;    % to store u in x direction, of w-th particle with time t 
            u(2,1,w,t)=u(2,1,w,t-1)+a(2,1,t)*t;
            u(3,1,w,t)=u(3,1,w,t-1)+a(3,1,t)*t;
            
%             end
    
            
%             fprintf('%u\n',u);

            s(1,1,w,t)=u(1,1,w,t)*t+(((1/2)*a(1,1,w))*t^2);   % s=ut+1/2at^2;
            s(2,1,w,t)=u(2,1,w,t)*t+(((1/2)*a(2,1,w))*t^2);
            s(3,1,w,t)=u(3,1,w,t)*t+(((1/2)*a(3,1,w))*t^2);
            
            
%             ds=s_fin-s_ini;
            ds(1,1,w,t)=s(1,1,w,t)-s(1,1,w,t-1);
            ds(2,1,w,t)=s(2,1,w,t)-s(2,1,w,t-1);
            ds(3,1,w,t)=s(3,1,w,t)-s(3,1,w,t-1);

            x_change(w,t)=ds(1,1,w,t);
            y_change(w,t)=ds(2,1,w,t);
            z_change(w,t)=ds(3,1,w,t);
            
            x_final(w,1)=0;
            y_final(w,1)=0;
            z_final(w,1)=0;
            
            x_final(w,t)=x_final(w,t-1)+x_change(w,t);

            y_final(w,t)=y_final(w,t-1)+y_change(w,t);

            z_final(w,t)=z_final(w,t-1)+z_change(w,t);
            
            x(w)=x_final(w,t);
            y(w)=y_final(w,t);
            z(w)=z_final(w,t);
            
%             fprintf('%z\n',z(w));
        
% %             
%             scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
%             pause(0.1);
            
            
            
            
%             NEED TO UPDATE particledata / 3D to change accordingly
%             HOW TO?
%             ballp(i)=(x(i),y(i),z(i))
%             dt = 0.5
%             t=0.0
% 
% 
%             while True:
% 
%               t = t + dt
%               ballpos(i) = ballpos(i) + (ballp(i)/m(i))*dt
% 
%               if not (side > ball.x > -side):   (0,Y,Z) OR (100,Y,Z)
%                 ball.p.x = -ball.p.x
%               if not (side > ball.y > -side):   (X,0,Z) OR (X,100,Z)
%                 ball.p.y = -ball.p.y
%               if not (side > ball.z > -side):
%                 ball.p.z = -ball.p.z   HERE BALL IS STUCK TO TOP N BOTTOM  for any (X,Y,0) or (X,Y,100)

  

    end
            
end

% 
% for w=1:1000
%     for t=2:1:10
% 
%     scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
%     pause(0.1);
%     end
% end
%  
%

 

[ApexOfDE]

I need help updating these new x,y,z postion at time range (1:10), I am sort of stuck there.

I don't understant what you means :( sry. When I glance at your code, I see that you can find x/y/z_final. You just need to draw a figure again to obtain new pos.

I also need to lock down this 3D space to a 100x100m box, it seems to go beserk,...

do you mean you want to see fig in 2d? If that, view(2) can help you.

these particles should get stuck to the top and bottom if they reach there... and if it goes beyond left or right, it should appear in the opposite side.

It calls periodic condition. You can do this simply. Let's assume that z-pos ranges from 0 to 100 and there is another matrix to indicate whether particles have its z-pos out of this range (0 - no, 1 - yes). If (1), skip calculating this particle.

For left-right side, you can put a condition after calculating and modifying data for this special move. It should be more complicated.

p/s: I can only say this because I don't have time to read your code and completely understand it. Ill try it again when I am free :D

 

[logix88]

I don't understant what you means :( sry. When I glance at your code, I see that you can find x/y/z_final. You just need to draw a figure again to obtain new pos.



do you mean you want to see fig in 2d? If that, view(2) can help you.



It calls periodic condition. You can do this simply. Let's assume that z-pos ranges from 0 to 100 and there is another matrix to indicate whether particles have its z-pos out of this range (0 - no, 1 - yes). If (1), skip calculating this particle.

For left-right side, you can put a condition after calculating and modifying data for this special move. It should be more complicated.

p/s: I can only say this because I don't have time to read your code and completely understand it. Ill try it again when I am free :D

Let me elaborate... I have x,y,z position of a 1000 particles. Now, over time these particles will be moving, due to the different forces applied on them, I want to simulate how these particles move around this 100x100x100m box. The t(1:10) was just time range going form 1 to 10, in steps of 1.

the 100mx100m was a mistake, it is supposed to be 100x100x100m

I have updated the code, to the one previously used, to make things a little clear.

Thanx!

 

[ApexOfDE]

It is a simple code that i tried to write. In this code, I used your coordinate data and created new position by using random integers. If z-pos is out of [0 100] range, this particle will not move next time.

clc
clear all

%
data = load('particledata.dat');
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
a_z = zeros(length(z), 1);

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

for t = 1 : 100
    for j = 1 : length(x)
        if a_z(j) == 1
            continue
        end
        
        d = randi(10, 1, 3);
        s = randi(2, 1, 3);
        
        x(j) = x(j) + d(1) * (-1)^s(1);
        y(j) = y(j) + d(2) * (-1)^s(2);
        z(j) = z(j) + d(3) * (-1)^s(3);
        
        if x(j) <= 0 || x(j) >= 100
            C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
            x(j) = 100 * C;
        end
        
        if y(j) <= 0 || y(j) >= 100
            C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
            y(j) = 100 * C;
        end
        
        %if p is at bottom or top, it can't move anymore.
        if z(j) < 0 || z(j) > 100
            C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
            z(j) = 100 * C;
            a_z(j) = 1;
        end
    end
    
    scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
    pause(1)
end

Hope this helps.

 

[logix88]

It is a simple code that i tried to write. In this code, I used your coordinate data and created new position by using random integers. If z-pos is out of [0 100] range, this particle will not move next time.

clc
clear all

%
data = load('particledata.dat');
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
a_z = zeros(length(z), 1);

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

for t = 1 : 100
    for j = 1 : length(x)
        if a_z(j) == 1
            continue
        end
        
        d = randi(10, 1, 3);
        s = randi(2, 1, 3);
        
        x(j) = x(j) + d(1) * (-1)^s(1);
        y(j) = y(j) + d(2) * (-1)^s(2);
        z(j) = z(j) + d(3) * (-1)^s(3);
        
        if x(j) <= 0 || x(j) >= 100
            C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
            x(j) = 100 * C;
        end
        
        if y(j) <= 0 || y(j) >= 100
            C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
            y(j) = 100 * C;
        end
        
        %if p is at bottom or top, it can't move anymore.
        if z(j) < 0 || z(j) > 100
            C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
            z(j) = 100 * C;
            a_z(j) = 1;
        end
    end
    
    scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
    pause(1)
end

Hope this helps.

Thnx for the code, Its seems pretty straight forward with random x,y,z positions, its just when trying to implement forces and mapping them wrt to time, it seemed a bit too confusing... I will give it a ago, and come bak if I have any problems, Thanx a lot for ur help!

 

[D H]

logix88,

Your code has a *big* problem in the way you are propagating the particles from one time step to the next. If you want realistic physics, those particles need to have position and velocity. The acceleration changes the velocity, and the velocity in turn changes the position.

Make sure you do things in that order. If you update the position first and velocity second you will have what is called the Euler integration technique. Simply switch the order (velocity, then position) and you will have what is called symplectic Euler integration technique. The difference between the two is immense. Basic Euler integration is truly abysmal. Symplectic Euler, while not that great, represents a vast improvement in realism.

 

[logix88]

logix88,

Your code has a *big* problem in the way you are propagating the particles from one time step to the next. If you want realistic physics, those particles need to have position and velocity. The acceleration changes the velocity, and the velocity in turn changes the position.

Make sure you do things in that order. If you update the position first and velocity second you will have what is called the Euler integration technique. Simply switch the order (velocity, then position) and you will have what is called symplectic Euler integration technique. The difference between the two is immense. Basic Euler integration is truly abysmal. Symplectic Euler, while not that great, represents a vast improvement in realism.

I thought that's what I did, first I calculate Forces due to heat, wind, gravity for each particle, and I get acceleration from those values. Initial vel is 0, I then use the acceleration to obtain u=u0+at and then plug that into s=ut+1/2*at^2. Then obtain the change in distance for each x,y,z between previous time step and current time step and finally update the particles position. Do you mean, that I should update final position before distance? What other way can I do that?

clc
clear all

%
data = load('particledata.dat');
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
m = data(:, 4);
a_z = zeros(length(z), 1);
g = 9.81;
q = 0.08;

A=data(:,5);

F_w=zeros(3,1,1000);

F_g=zeros(3,1,1000);

F_h=zeros(3,1,1000);

u=zeros(3,1,1000,100);

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

for t = 2 : 100
    for j = 1 : length(x)
%         if a_z(j) == 1
%             continue
%         end
        
        
        F_g_z(j)=-m(j)*g*z(j);

        F_g(1,1,j)=0;
        F_g(2,1,j)=0;
        F_g(3,1,j)=F_g_z(j);

        rx=50-y(j);
        ry=x(j)-50;

        F_w_x(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*rx;
        F_w_y(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*ry;
        F_w(1,1,j)=F_w_x(j);
        F_w(2,1,j)=F_w_y(j);
        F_w(3,1,j)=0;

        T(j)=1/16*(70*x(j)+(50-y(j))^2);

        F_h1(j)=q*A(j)*T(j)*z(j);

        F_h(1,1,j)=0;
        F_h(2,1,j)=0;
        F_h(3,1,j)=F_h1(j);

        F=F_g+F_w+F_h;
    %     
    %     a=a_g+a_w+a_h;
        a=F/m(j);

        u(1,1,j,1)=a(1,1,j);
        u(2,1,j,1)=a(2,1,j);
        u(3,1,j,1)=a(3,1,j);

        u(1,1,j,t)=u(1,1,j,t-1)+a(1,1,j)*t;
        u(2,1,j,t)=u(2,1,j,t-1)+a(2,1,j)*t;
        u(3,1,j,t)=u(3,1,j,t-1)+a(3,1,j)*t;

%             end


%             fprintf('%u\n',u);

        s(1,1,j,t)=u(1,1,j,t)*t+(((1/2)*a(1,1,j))*t^2);
        s(2,1,j,t)=u(2,1,j,t)*t+(((1/2)*a(2,1,j))*t^2);
        s(3,1,j,t)=u(3,1,j,t)*t+(((1/2)*a(3,1,j))*t^2);


%             ds=s_fin-s_ini;
        ds(1,1,j,t)=s(1,1,j,t)-s(1,1,j,t-1);
        ds(2,1,j,t)=s(2,1,j,t)-s(2,1,j,t-1);
        ds(3,1,j,t)=s(3,1,j,t)-s(3,1,j,t-1);

        x_change(j,t)=ds(1,1,j,t);
        y_change(j,t)=ds(2,1,j,t);
        z_change(j,t)=ds(3,1,j,t);

        x_final(j,1)=0;
        y_final(j,1)=0;
        z_final(j,1)=0;

        x_final(j,t)=x_final(j,t-1)+x_change(j,t);

        y_final(j,t)=y_final(j,t-1)+y_change(j,t);

        z_final(j,t)=z_final(j,t-1)+z_change(j,t);

        x(j)=x_final(j,t);
        y(j)=y_final(j,t);
        z(j)=z_final(j,t);
        
                
        if x(j) <= 0 || x(j) >= 100
            C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
            x(j) = 100 * C;
        end
        
        if y(j) <= 0 || y(j) >= 100
            C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
            y(j) = 100 * C;
        end
        
        %if p is at bottom or top, it can't move anymore.
        if z(j) < 0 || z(j) > 100
            C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
            z(j) = 100 * C;
            a_z(j) = 1;
        end
    end
    
    scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
    pause(0.1)
end

 

[D H]

I thought that's what I did, first I calculate Forces due to heat, wind, gravity for each particle, and I get acceleration from those values. Initial vel is 0, I then use the acceleration to obtain u=u0+at and then plug that into s=ut+1/2*at^2.

This assumes constant acceleration. You don't have constant acceleration.

What you should be doing is something like this:

\aligned<br /> \vec u(t+\Delta t) &amp;= \vec u(t) + \vec a(t)\Delta t \\<br /> \vec x(t+\Delta t) &amp;= \vec x(t) + \vec u(t+\Delta t) \Delta t<br /> \endaligned

You going from t=2 to 100 (and why 2?) in steps of 1. So in your case, \Delta t = 1. Do not do this:

\aligned<br /> \vec u(t+\Delta t) &amp;= \vec u(t) + \vec a(t)\Delta t \\<br /> \vec x(t+\Delta t) &amp;= \vec x(t) + \vec u(t)\Delta t<br /> \endaligned

That second one is the basic Euler method. The first one is symplectic Euler.

 

[logix88]

This assumes constant acceleration. You don't have constant acceleration.

What you should be doing is something like this:

\aligned<br /> \vec u(t+\Delta t) &amp;= \vec u(t) + \vec a(t)\Delta t \\<br /> \vec x(t+\Delta t) &amp;= \vec x(t) + \vec u(t+\Delta t) \Delta t<br /> \endaligned

You going from t=2 to 100 (and why 2?) in steps of 1. So in your case, \Delta t = 1. Do not do this:

\aligned<br /> \vec u(t+\Delta t) &amp;= \vec u(t) + \vec a(t)\Delta t \\<br /> \vec x(t+\Delta t) &amp;= \vec x(t) + \vec u(t)\Delta t<br /> \endaligned

That second one is the basic Euler method. The first one is symplectic Euler.

u(1,1,j,1)=a(1,1,j);
u(2,1,j,1)=a(2,1,j);
u(3,1,j,1)=a(3,1,j);

u(1,1,j,t+1)=u(1,1,j,t)+a(1,1,j)*t;
u(2,1,j,t+1)=u(2,1,j,t)+a(2,1,j)*t;
u(3,1,j,t+1)=u(3,1,j,t)+a(3,1,j)*t; x_final(j,1)=x(j);
y_final(j,1)=y(j);
z_final(j,1)=z(j);

x_final(j,t+1)=x_final(j,t)+u(1,1,j,t+1);
y_final(j,t+1)=y_final(j,t)+u(2,1,j,t+1);
z_final(j,t+1)=z_final(j,t)+u(3,1,j,t+1);

x(j)=x_final(j,t+1);
y(j)=y_final(j,t+1);
z(j)=z_final(j,t+1);

Is that what it should be like?

 

[D H]

Is that what it should be like?

No. Why are you doing this?

u(1,1,j,t+1)=u(1,1,j,t)+a(1,1,j)*t;
u(2,1,j,t+1)=u(2,1,j,t)+a(2,1,j)*t;
u(3,1,j,t+1)=u(3,1,j,t)+a(3,1,j)*t;

Think of it this way: Suppose you had time going from 102 to 200? None of your forces depend on time. When time is zero is just a label. Moving that label doesn't change the physics.

 

[logix88]

No. Why are you doing this?

Think of it this way: Suppose you had time going from 102 to 200? None of your forces depend on time. When time is zero is just a label. Moving that label doesn't change the physics.

Hey, I am sry I don't understand what you said, I am not a physicist, I knw I have some prb understanding the simple physics involved in this.

TBH, I am more confused now than before...