Charges in constant electromagnetic fields

This code calculate the trajectory and the velocity of a charge in a constant electric and magnetic field.

The problem is the following:

m\frac{d\vec{v}}{dt} = q(\vec{E} + \vec{v}\times\vec{B})
\frac{d\vec{x}}{dt} = \vec{v}

Where mass, charge and velocity are in non-dimensional form:

m = M/M_{electron}, q = Q/e, v = V/c

The easiest way to solve this problem is by using euler forward method, of course it have several problems:

— it’s first order

— can be unstable

— energy is not conserved

So maybe it’s better to use centered finite differences that is second a order method and conserve energy. In order to do that we need to use two different time discratization for velocity and position (leapfrog). The discretized equation become:

m\frac{\vec{v}^{n+1/2} - \vec{v}^{n-1/2} }{\Delta{t}} = q(\vec{E} + \vec{v}^{n}\times\vec{B})
\frac{\vec{x}^{n+1} - \vec{x}^n}{\Delta{t}} = \vec{v}^{n+1}

Where whe can approximate:

\vec{v}^n = \frac{\vec{v}^{n+1/2}+\vec{v}^{n-1/2}}{\Delta t}

To solve this we must first put to the left the terms in (n+1/2) and to the right the terms in (n-1/2). Then we can build the matrices needed to solve the system of equations. (in the code the matrices are called A and M)

The following is the Matlab/Octave code:

close all
clear
clc
%LEAPFROG METHOD
%v(n+1/2) - v(n+1/2)/2c x B = v(n-1/2) + 1/m(E + v(n-1/2)/2c x B )

%particle characteristics - nondimensional
q = -1;						%charge
m = 1;						%mass
 
%time definition
t0 = 0;
tf = 1;
Nt = 100;
dt = (tf-t0)/Nt;
tx = linspace(t0,tf,Nt);  	%space time discretization
tv = tx + dt/2; 			%velocity  time discretization 
%plot(tx,0,'o')

%MATRIX DEFINITION

V = zeros(3,Nt);	%velocity matrix
X = zeros(3,Nt);	%coordinate matrix

%IC

V(1,1) = 0;
V(2,1) = 1;
V(3,1) = 0;
X(1,1) = 0;
X(2,1) = 0;
X(3,1) = 0;

%CONSTANT FIELD DEFINITION

B = zeros(3,1);
E = zeros(3,1);

B(1) = 10;
B(2) = 0;
B(3) = 0;

E(1) = 1;
E(2) = 0;
E(3) = 0;

%SYSTEM MATRIX
%Av(n+1/2) = Mv(n) + dt*q/m*E
A = [1, -dt*q/m*B(3), + dt*q/m*B(2);
	 dt*q/m*B(3), 1, -dt*q/m*B(1);
	 -dt*q/m*B(2), dt*q/m*B(1), 1];
	 
	 
M = [1, dt*q/m*B(3), -dt*q/m*B(2);
	 -dt*q/m*B(3), 1, +dt*q/m*B(1);
	 dt*q/m*B(2), -dt*q/m*B(1), 1];	 
	 
	 
%SOLUTION CYCLE	 
for n = 1:Nt-1
	
	V(:,n+1) = A\(M*V(:,n) + dt*q/m*E);
		
	%VMOD = norm(V(:,n+1));
		
	X(:,n+1) = X(:,n) + dt*V(:,n+1);
		
	%XMOD = norm(X(:,n+1));
		
end

%PLOT

plot3(X(1,:),X(2,:),X(3,:))
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')
	 

• 

•