MatLab: Motion of charged particle in electromagnetic field in cylinder

Question

We wish to solve the differential equations for a particle's movement in a electromagnetic field inside a cylinder. This has to be done using the a non built-in Runge-Kutta method of the 4th order. We are given the initial velocity v0, magnetic constant B, electric constant E, the cylinder radius R, the midpoint (a,b) of the cylinder, the weight $\omega$, the charge Q, and the particle mass m. The magnetic field is B = B*[0, 0, 1 - $\omega \frac{r_1}{R}$], where $r1 = \sqrt{ (x - a)^2 + (y - b)^2 }$, and is pointed in the z-direction. The electric field is given by E= E*[1, 0, 0]. Given that the velocity of the particle is 0 in the z-direction outside the cylinder, the movement of the particle in the cylinder is given by the 2nd order differential equations $ \ddot x = \frac{Q}{m}(E + B\dot y(1 - w \frac{r_1}{R})) $ and $ \ddot y = -\frac{Q}{m}B\dot x(1 - w \frac{r_1}{R}) $, which have been derived from $ m \ddot r = Q($ E $ + \ \dot r \times $ B$)$.

$\ddot x$ and $\ddot y$ can be reduced to a system of 1st order differential equations ...

function dudt = uprim(t, u)
    m = 9.1091*10^-31; 
    Q = -1.6021*10^-19;  
    E = 20.0;            
    B = 0.92*10^-4;      
    R = 0.012;           
    a = 0.015;         
    b = 0;             
    w = 0.2;             

    x = u(1);
    xDot = u(2);

    y = u(4);
    yDot = u(5);

    xDotDot =  (Q/m)*(E + B*u(5)*(1 - w*(sqrt( (u(1) - a).^2 + (u(4) - b).^2 )/R)));
    yDotDot =      -(Q/m)*B*u(2)*(1 - w*(sqrt( (u(1) - a).^2 + (u(4) - b).^2 )/R));

    dudt = zeros(size(u));
    dudt(1) = xDot;
    dudt(2) = xDotDot;
    dudt(3) = yDot;
    dudt(4) = yDotDot;
end

The particle moves linearly in the positive direction along the x-axis, and enters the cyliner at (a - R, 0). Using initial values and settings for the Runge-Kutta method we get ...

v0 = 455*10^3;
h = 0.001; 
t = 0:h:1;
y = zeros(2, length(t));`
y(1, 1) = a - R; 
y(2, 1) = v0; 
y(3, 1) = 0; 
y(4, 1) = 0;
y(5, 1) = 0;     
y(6, 1) = 0;
for i = 1:(length(t) - 1)
    k1 = uprim(t(i)        , y(:, i)           );
    k2 = uprim(t(i) + 0.5*h, y(:, i) + 0.5*h*k1);
    k3 = uprim(t(i) + 0.5*h, y(:, i) + 0.5*h*k2);   
    k4 = uprim(t(i) +     h, y(:, i) +     h*k3);
    y(:, i + 1) = y(:, i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h;
end

We're stuck here. The values of y are either NaN or infinitely high or low. Maybe we could use some form or parametization for x and y, but were aren't getting any further on that either. Can we get some hints on what we should do, or some pointers to where the code is incorrect? Thanks

Answer 1

You have to take the scale of your constants into account. The Lipschitz constant L of the ODE system can be taken as abs(Q/m*B)=1.6180874e+07. For the RK4 method to work sensibly you need L*h between 1e-4 and 1e-2, for the larger choice this is h=1e-9. The upper bound has then to be adapted to be conform with the number of integration steps, below I used 10000 steps to see a significant portion of the solution curve.

Next correct the use of indices, if yDot and yDotDot are on position 3 and 4 in the output dudt of uprime, then y and yDot have also to be found at these positions in the state vector.

R = 0.012;           
a = 0.015;         
b = 0;             
w = 0.2;     
v0 = 455*10^3;

function dudt = uprim(t, u)
    m = 9.1091*10^-31; 
    Q = -1.6021*10^-19;  
    E = 20.0;            
    B = 0.92*10^-4;      

    x = u(1);
    xDot = u(2);

    y = u(3);
    yDot = u(4);

    r = sqrt( (x - a).^2 + (y - b).^2 );
    xDotDot =  (Q/m)*(E + B*yDot*(1 - w*(r/R)));
    yDotDot =  (Q/m)*(  - B*xDot*(1 - w*(r/R)));

    dudt = zeros(size(u));
    dudt(1) = xDot;
    dudt(2) = xDotDot;
    dudt(3) = yDot;
    dudt(4) = yDotDot;
end

h = 1e-9; 
t = 0:h:1e-5;
u = zeros(6, length(t));
u(1, 1) = a - R; 
u(2, 1) = v0; 
u(3, 1) = 0; 
u(4, 1) = 0;
u(5, 1) = 0;     
u(6, 1) = 0;
for i = 1:(length(t) - 1)
    k1 = uprim(t(i)        , u(:, i)           );
    k2 = uprim(t(i) + 0.5*h, u(:, i) + 0.5*h*k1);
    k3 = uprim(t(i) + 0.5*h, u(:, i) + 0.5*h*k2);   
    k4 = uprim(t(i) +     h, u(:, i) +     h*k3);
    u(:, i + 1) = u(:, i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

figure(1);
plot(t, u(1,:), t, u(3,:));
figure(2);
plot(u(1,:),u(3,:));

The plot contains left the x and y components as functions of t and right the curve (x,y).