Simulation of charged particle in matlab

35 次查看（过去 30 天）
simar
simar 2011-6-9
评论： Walter Roberson 2019-1-21
I made this code but it does not give correct result probably due to quantisation error .. any suggestions how can i fix it??
I want to have a simulation if this kind http://www.youtube.com/watch?v=a2_wUDBl-g8
% script that simulates a moving particle with some initial velocity in a
% magnetic field B
v0 = [5 0 0]; %initial velocity
B = [0 0 -5]; %magnitude of B
m = 5; % mass
q = 1; % charge on particle
r0 = [0 0 0]; % initial position of particle
t = 0;
% Now we want to find the next velocity as the particle enters the magnetic
% field and hence its new position
r = r0;
v = v0;
dt = 0.00000000000000000001;
figure
xlim([-25 25])
ylim([-25 25])
hold on
for n = 1:100
%plot it
plot(r(1),r(2),'*');
%pause
% update time
t = t+dt;
% new position r
dr = v*dt;
r = r + dr;
%find new velocity
dv = (q/m) * cross(v,B);
v = v + dv;
end
  3 个评论
显示 1更早的评论
Jan
Jan 2011-6-10
"Does not give correct results" is not a useful description of the problem. How could we know, what the correct result is to find the error?!
Anyhow, "dt = 0.00000000000000000001" is ridiculous small. Be aware than 1+dt==1!
simar
simar 2011-6-10
Andrew and Jan
I will note this in future ... Thanks

请先登录，再进行评论。

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

采纳的回答

Andrew Newell
Andrew Newell 2011-6-10
There is nothing wrong with the physics. I think it is just the numerical stability of your code. The MATLAB ODE suite handles this much better:
v0 = [5 0 0.1]'; %initial velocity
B = [0 0 -5]'; %magnitude of B
m = 5; mass
q = 1; charge on particle
r0 = [0 -5 0]'; initial position of particle
tspan = [0 100];
The code below solves the system of equations
x' = vx
y' = vy
z' = vz
vx' = (q/m) (v x B)_x
vy' = (q/m) (v x B)_y
vz' = (q/m) (v x B)_z
The approach is to treat the position and velocity as independent variables.
y0 = [r0; v0];
f = @(t,y) [y(4:6); (q/m)*cross(y(4:6),B)];
[t,y] = ode45(f,tspan,y0);
plot3(y(:,1),y(:,2),y(:,3))
  9 个评论
显示 7更早的评论
Walter Roberson
Walter Roberson 2019-1-16
The second parameter to f, namely y, will be a column vector, so y(4:6) will be a column vector.
cross(y(4:6), B) will be a column vector.
E is a row vector. Row vector plus column vector was an error before R2016b, but as of R2016b started being treated similarly to bsxfun. So your code is like
f = @(t,y) [y(4:6); (q/m)*( bsxfun(@plus, E, cross(y(4:6),B) ))];
The 1 x 3 row vector plus 3 x 1 column vector will produce a 3 x 3 result.
You are then trying to use [;] between a 3 x 1 from y, and the 3 x 3 from the remainder of the computation. That is going to fail because of the mismatch on the number of columns.
Solution: change your definition to
E = [0 15000 0].' ;

请先登录，再进行评论。

更多回答（6 个）

Ivan van der Kroon
Ivan van der Kroon 2011-6-10
You probably want to have something like
fun=@(t,x) [x(4:6);cross(x(4:6),q/m*B)];
t=linspace(0,tend,1e3);
[t,x]=ode45(fun,t,[r0;v0]);
For the movie read about getframe:
doc getframe
simar
simar 2011-6-10
Thanks a lot friends for your answers. I'm just a matlab beginner at present and is also learning mathematics.
I thnik in my method I used some approximation to solve differential equations. In some sense that taylor method, I think I realized that those are differential equations but I dinnt knew the way to implement that. Nevertheless.. thanks for your help
Here is how I coded. Pleas let me know any suggestions...
% script that simulates a moving particle with some initial velocity in a
% magnetic field B
v = [3 4 1]; %initial velocity
B = [0 0 -5]; %magnitude of B
m = 5; % mass
q = 1; % charge on particle
r0 = [5 0 0]; % initial position of particle
t = 0;
%find velocity parallel to B and perpendicular to B
v_para = (dot(v,B)/norm(B))*(B/norm(B));
v_per = v-v_para;
% find the circles centre
r = m*(norm(v_per))/(q*norm(B));
theta = atan(v_per(2)/v_per(1))+pi/2;
xc=r0(1)+r*cos(theta);
yc=r0(2)+r*sin(theta);
w= norm(v_per)/r;
figure
plot3(-10:0.1:10,0,0);
hold on;
plot3(0,-10:0.1:10,0);
plot3(0,0,-10:0.1:10)
xlim([-15 15])
ylim([-15 15])
%zlim([-15 15])
t=0;
%dt=0.01;
tic
for n=1:4000
dt = toc;
tic
x=xc+r*cos(w.*t + pi+theta);
y=yc+r*sin(w.*t + pi+theta);
z=v_para*t;
plot3(x,y,z,'--.');
hold on
t=t+dt;
pause(0.00000000005)
end
I hope I did it better this time
  4 个评论
显示 2更早的评论
simar
simar 2011-6-11
I was just thinking to add some transformatin code so that when I change B it can chang its axis .. of motion ..
can you help me with some source or code .....????
Andrew Newell
Andrew Newell 2011-6-11
Time to start a new question, I think.

请先登录，再进行评论。

Matt Tearle
Matt Tearle 2011-6-10
Without knowing the math behind what you're trying to simulate, I don't know if it's a problem in the way you've coded up the equations or if it's an issue of numerical implementation. But either way the culprit is in the line
dv = (q/m)*cross(v,B);
If you increase the number of loops or increase dt, you'll see that you get an outward spiral of points. If you keep a track of norm(dv), you'll see that this is because the magnitude of dv grows exponentially. So with even an absurdly small dt like you have, dr will eventually explode (and r with it). As a result, you'll see nothing for a while (a single point at the origin), then a very quick spiral outwards, before norm(r) just blows up to Inf.
IOW, you have a classic runtime error: MATLAB is doing exactly what you asked it to do. Just not, it seems, what you want it to do :)
  5 个评论
显示 3更早的评论
Andrew Newell
Andrew Newell 2011-6-10
The I Ching says it best:
Heaven is full of electrons
Spiralling round field lines.
Thus, the superior man studies
Maxwell's equations.

请先登录，再进行评论。

simar
simar 2011-6-16
Hey guys check out this one .. http://www.mathworks.com/matlabcentral/fileexchange/2268-projects-for-scientists-and-engineers/content/penland/lorentz.m Some has implemented the same result as mine and isn't getting the correct result. But how does then matlab able to calculate the correct path using its ODE solver. Which algorithm it applies. I tested it for large time still there is no error in trajectory...
  2 个评论
Bjorn Gustavsson
Bjorn Gustavsson 2011-6-16
Calculation of particle trajectories in magnetic fields is a tricky problem. Even for your case with a constant B all but one of the ODE-solvers of matlab fails. I ran Andrew Newell's example but multiplied the end time with 100, then you'll see that there is variation of the gyro radius with time - and it shouldn't be.
One problem with particle motions in magnetic fields is that the particle energy should be conserved - for your case it is as simple as conservation of the kinetic energy:
(lets use "'" for time derivatives)
v' = q/m*cross(v,B) ->
dot(v,v') = q/m*dot(v,cross(v,B)) = 0,
since cross(v,B) is perpendiculat to B.
integrate once with respect to time and you get:
v^2/2 = Const.
The only ODE-solver that conserves that is ode23t, the others don't, some increase v^2 some decrease. Search for Störmer or Verlet integration for more information.
Andrew Newell
Andrew Newell 2011-6-16
The documentation for the ODE suite (http://www.mathworks.com/help/techdoc/ref/ode23.html) identifies the algorithm for each function.

请先登录，再进行评论。

Suraj Tamrakar
Suraj Tamrakar 2017-7-20
编辑：Suraj Tamrakar 2017-7-20
Could someone help me clarify the math for the parallel and perp velocity component .
v_para = (dot(v,B)/norm(B))*(B/norm(B));
I don't quite get this formula used above. This is quite unfamiliar expression to me. Thank you
Mathew Orman
Mathew Orman 2019-1-21
You have wrong physics and your simulation should show gain of kinectic energy while charge is accelerated by electric force field...
https://www.youtube.com/watch?v=lhldn0ef138&feature=youtu.be
  1 个评论
Walter Roberson
Walter Roberson 2019-1-21
Samir comments:
The reason for the increase in kinetic energy is numerical errors in the simulation. By the way, dynamo of my bicycle works just fine.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by