Leapfrog integration with multiple particles

Hi all i'm new to programming and I wanted to ask a quick question: Consider the following leapfrog integration for one particle:

time=100;

steps=1000000;

dt=time/(steps-1);

v0=5

theta=2*pi*rand(1);

vx(1)=v0*sin(theta);

vy(1)=v0*cos(theta);

t(1)=0;

x(1)=50;

y(1)=50;

for i=1:steps-1

      t(i+1)=t(i)+dt;
      x(i+1)=x(i)+vx(1)*dt;
      y(i+1)=y(i)+vy(1)*dt;

end

How do I increase the number of particles in the leapfrog integration? eg 1000? Thanks in advance!

Some info (edited after request): Ok so, basically the particle starts at time t=0, at position x,y and at velocity vx,vy.

The simulation lasts for 100 seconds and consists in a million step. Leapfrog integration, in this case, consists in the for loop that updates the particle position every 100/1000000.

The problem is that i dont know how to scale this for loop for multiple particles. Thanks!

[Geoff Hayes]

Frederico - how is the single particle used in the above code? I guess some background on leapfrog integration would help. :)

[AFHG]

Ok so, basically the particle starts at time t=0, at position x,y and at velocity vx,vy.

The simulation lasts for 100 seconds and consists in a million step. Leapfrog integration, in this case, consists in the for loop that updates the particle position every 100/1000000.

The problem is that i dont know how to scale this for loop for multiple particles. Thanks!

Answer 1

[Star Strider]

This works (with a much smaller value for ‘steps’):

time=100;

steps=1000000;
% steps = 100;

dt=time/(steps-1);

v0=5;

NrParticles = 1000;

theta=2*pi*rand(1,NrParticles);

vx=v0*sin(theta);

vy=v0*cos(theta);

t(1)=0;

x(1,1:NrParticles)=50;

y(1,1:NrParticles)=50;

for k = 1:NrParticles
    for i=1:steps-1

          t(i+1)=t(i)+dt;
          x(i+1,k+1)=x(i,k)+vx(i)*dt;
          y(i+1,k+1)=y(i,k)+vy(i)*dt;

      end
  end

figure
hold all
for k1 = 1:NrParticles
    plot(x(:,k1), y(:,k1))
end
hold off

You will likely have to experiment with your code. I changed it slightly to use different values for ‘vx’ and ‘vy’ in the loop.

EDIT — Corrected error in ‘k’ loop.

[Star Strider]

I already added that in my original Answer. It is the ‘k1’-loop after the ‘k’ loop. (I usually plot the results of my Answers, where appropriate, to check my code.)

The ‘k’-loop edit did not affect the plot call.

[AFHG]

Try to run the code it just comes weird on the plot.

[AFHG]

and, if you look into the y values that are updated for the same particles they become negative, positive, negative etc, which is impossible for constant velocity. This happens the second time the values of x and y are updated.

[Star Strider]

I provided a solution to scaling up your code for multiple particles. That is what you asked.

As I understand your code, the x- and y-values are position coordinates, so they could wander randomly over the space you define for them, both positive and negative, constant velocity or not.

The rest I defer to you to troubleshoot, since I have no idea what you are doing.

[AFHG]

Yes, I know that but if you look into them they go from 50, 47, 0.7,0.4,etc and all stay around 0. This means that the leapfrog integration is updating wrongly the y and x values and was thus not scaled appropriately for more particles. As you prob know much more about matlab than me since i started using it yesterday I kindly ask for your help to fix this. Thanks

[Star Strider]

Your code (and my additions) run without error.

I know nothing about Leapfrog integration (link) other than what is in this article. However, you defined velocity (that I assume are ‘vx’ and ‘vy’) to have varying signs, since you defined them using trigonometric functions of a complete circle:

theta=2*pi*rand(1,NrParticles);
vx=v0*sin(theta);
vy=v0*cos(theta);

If you want the velocities to be only positive, you have to define them that way.

[AFHG]

Negative velocities are not a problem, my point is: if a particle starts at 50 and has a negative velocity of -3, its position in time will be: 50, 47, 44... and eventually become negative. the same particle will never become positive (again) or stationary around 0 as the vector is constant. It is also impossible for a velocity of less than 5 to go from 47 to 0 in one time-step. There must be a mistake in the loop. I attached as a png, a screen of the values i get

[Star Strider]

You must troubleshoot that. All I did was to scale up your code for multiple particles. I did not otherwise change your original code.

[AFHG]

My original code worked this doesn't but ok thank you anyways i will try find a solution

[AFHG]

I just managed to fix it thank you for everything

[AFHG]

This was the correct loop, btw:

e
for k = 1:NrParticles
  for i=1:steps-1
        t(i+1)=t(i)+dt;
        x(i+1,k)=x(i,k)+vx(1,k)*dt;
        y(i+1,k)=y(i,k)+vy(1,k)*dt;
    end
end

[Star Strider]

As always, my pleasure.

I was not certain how to scale up ‘vx’ and ‘vy’.