Simulation of charged particle in matlab

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

[Andrew Newell]

Please read http://www.mathworks.com/matlabcentral/answers/7885-tutorial-how-to-format-your-question.

[Jan]

"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]

Andrew and Jan

I will note this in future ... Thanks

Answer 1

[Andrew Newell]

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))

[simar]

It works thanks ..

[simar]

I figured out its the numerical stability of the code. Its the quantisation error. I think.

BTW I think I applied taylor method (in some sense) to solve the DEs but the error was sgnificant because I was dealing with vectors. But how does ODE manage to so this without quantisation error??

[Bjorn Gustavsson]

For this equation ODE45 has some unwanted "numerical energy dissipation", ODE23t gives better results (conserves gyro-radius) if one integrates the equation over a factor 100 longer time span.

[simar]

ok i tried a lot, but I could not make out how can you define a function like this

f = @(t,y) [y(4:6); (q/m)*cross(y(4:6),B)];

I know this anonymous function but in that case only a function like

fun = @(x) (sin(x)) makes sense to me ..??

[Andrew Newell]

Have you read the documentation about anonymous functions? Like any function, it can accept a vector input and give a vector output. The expression [y(4:6); (q/m)*cross(y(4:6),B)] combines two vectors of length 3 to make a vector of length 6 (as long as y is a column vector).

Unlike a function in a separate file, it can treat variables like q, m and B as parameters; their values are fixed while y is allowed to vary. See the example at http://www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html#f4-70238.

I have tested the code in my answer and it works if you just cut and paste it into your command window.

[simar]

ok now I got it. Its how an anonymous function accepts vector inputs and vector outputs. But I have seen it before, even in documentation.

[Samir]

Well, i am whole 7 years late on this but could any one probably tell me how to implement an electric field E in the above mentioned script. I tried this but fails.

v0 = [5 0 0]';  %initial velocity
B = [0 0 0.5]';  %magnitude of B
E = [0 15000 0]
vd=E(2)/B(3)
v0 = [vd 0 0]';  %initial velocity
m = 40*1.6e-27;  %mass
q = 1.6e-19;  %charge on particle
r0 = [0 0 0]'; %initial position of particle 
tspan = [0:1e-6:1e-4];
y0 = [r0; v0];
f = @(t,y) [y(4:6); (q/m)*( E + cross(y(4:6),B) )];
[t,y] = ode45(f,tspan,y0);
plot3(y(:,1),y(:,2),y(:,3))
figure
plot(y(:,1),y(:,2))

Error with some dimentional missmatch.

[Walter Roberson]

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].' ;

[Samir]

Thanks, now i realize.

Answer 2

[Ivan van der Kroon]

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]

hmm ya getframe will help .. thanks a lot

Answer 3

[simar]

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

[Andrew Newell]

So you are plotting the analytical solution. That's good. It is strange, though, to use tic and toc, which are timing functions, to determine dt. Why not just use a fixed dt? Finally, pausing for 0.05 nanoseconds won't have any discernible effect.

[simar]

oh yeah... I wanted to make the particle move as if it were move if we cary an experiment .. but it moved so fast that i had to add a delay .. So tic toc can be removed though but i dinnt bothered cos its working fine ..

[simar]

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]

Time to start a new question, I think.

Answer 4

[Matt Tearle]

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 :)

[Matt Tearle]

Oooohhhh! Just read Andrew's answer. Yes, what he said. I (being dense) didn't realize you were doing Euler's method on ODEs. That's why you get the behavior I was describing. Do what Andrew says, and you're golden.

[Andrew Newell]

Thanks, Matt. FYI: (q/m)*cross(v,B) is the Lorenz force on a particle with mass m and charge q due to a magnetic field B.

[Matt Tearle]

Hmm, so I probably knew that for about 7.3 minutes. 17 years ago :)

[Andrew Newell]

The I Ching says it best:

Heaven is full of electrons

Spiralling round field lines.

Thus, the superior man studies

Maxwell's equations.

[simar]

Nice eplanation .. thanks a lot

Answer 5

[simar]

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...

[Bjorn Gustavsson]

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]

The documentation for the ODE suite (http://www.mathworks.com/help/techdoc/ref/ode23.html) identifies the algorithm for each function.

Answer 6

[Suraj Tamrakar]

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

Answer 7

[Mathew Orman]

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

[Walter Roberson]

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.