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

In summary, the code simulates a 3D particle simulation by considering forces acting on the particles (gravitation, wind, and heat). Updates to the x, y, z positions are made at time intervals (1:10). A force is applied to each particle to determine its acceleration. The code also updates the position of the camera (using the camera's position and the current time).
  • #1
logix88
12
0
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

Code:
%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
%  
%
 
Last edited by a moderator:
Physics news on Phys.org
  • #2
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
 
  • #3
ApexOfDE said:
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!
 
  • #4
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.

Code:
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.
 
  • #5
ApexOfDE said:
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.

Code:
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!
 
  • #6
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.
 
  • #7
D H said:
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?


Code:
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
 
  • #8
logix88 said:
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:

[tex]\aligned
\vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
\vec x(t+\Delta t) &= \vec x(t) + \vec u(t+\Delta t) \Delta t
\endaligned[/tex]

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

[tex]\aligned
\vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
\vec x(t+\Delta t) &= \vec x(t) + \vec u(t)\Delta t
\endaligned[/tex]

That second one is the basic Euler method. The first one is symplectic Euler.
 
  • #9
D H said:
This assumes constant acceleration. You don't have constant acceleration.

What you should be doing is something like this:

[tex]\aligned
\vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
\vec x(t+\Delta t) &= \vec x(t) + \vec u(t+\Delta t) \Delta t
\endaligned[/tex]

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

[tex]\aligned
\vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
\vec x(t+\Delta t) &= \vec x(t) + \vec u(t)\Delta t
\endaligned[/tex]

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?
 
Last edited:
  • #10
logix88 said:
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.
 
  • #11
D H said:
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...
 

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

1. What is particle simulation in 3D?

Particle simulation in 3D is a scientific technique used to model the behavior and movement of individual particles in a three-dimensional space. These particles can represent anything from atoms and molecules to stars and galaxies, and their interactions are simulated using mathematical equations and algorithms.

2. What is the purpose of particle simulation in 3D?

The purpose of particle simulation in 3D is to study and understand the complex behavior of systems with large numbers of individual particles. It allows scientists to predict how these particles will interact and behave under different conditions, providing valuable insights into various natural and artificial phenomena.

3. What are some common applications of particle simulation in 3D?

Particle simulation in 3D has a wide range of applications in fields such as physics, chemistry, biology, and engineering. Some common uses include studying the behavior of fluids, analyzing the dynamics of gases, simulating the formation and evolution of galaxies, and designing new materials with specific properties.

4. What tools and techniques are used for particle simulation in 3D?

Particle simulation in 3D typically involves using advanced computer programs and algorithms to simulate the interactions between particles. These programs may use techniques such as molecular dynamics, Monte Carlo simulations, or lattice-based methods to model the particles and their movements accurately.

5. What are the benefits of using particle simulation in 3D?

Particle simulation in 3D offers several benefits, including the ability to study complex systems that are difficult or impossible to observe in real life, the ability to simulate and test various scenarios quickly and efficiently, and the potential to discover new insights and make predictions about natural and artificial systems.

Similar threads

  • MATLAB, Maple, Mathematica, LaTeX
Replies
4
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
  • Programming and Computer Science
Replies
2
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
3K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
398
  • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
1K
Replies
2
Views
850
  • MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
2K
  • Advanced Physics Homework Help
Replies
1
Views
858
Back
Top