Leapfrog integration with multiple particles

2018 9 8
1 个回答
回答已采纳
更新时间：2018 9 8
4 次查看（30 天）

0 个投票

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!

2 个评论
显示 无 隐藏 无

Geoff Hayes
Geoff Hayes 2018-9-8
Frederico - how is the single particle used in the above code? I guess some background on leapfrog integration would help. :)
AFHG
AFHG 2018-9-8
编辑：AFHG 2018-9-8
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!

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Star Strider
Star Strider 2018-9-8
编辑：Star Strider 2018-9-8

0 个投票

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.

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Star Strider
Star Strider 2018-9-8
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
AFHG 2018-9-8
Try to run the code it just comes weird on the plot.
AFHG
AFHG 2018-9-8
编辑：AFHG 2018-9-8
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
Star Strider 2018-9-8
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
AFHG 2018-9-8
编辑：AFHG 2018-9-8
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
Star Strider 2018-9-8
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
AFHG 2018-9-8
编辑：AFHG 2018-9-8
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
Star Strider 2018-9-8
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
AFHG 2018-9-8
编辑：AFHG 2018-9-8
My original code worked this doesn't but ok thank you anyways i will try find a solution
AFHG
AFHG 2018-9-8
I just managed to fix it thank you for everything
AFHG
AFHG 2018-9-8
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
Star Strider 2018-9-8
As always, my pleasure.
I was not certain how to scale up ‘vx’ and ‘vy’.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

产品

版本

R2018a

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by