Una partícula de masa m y carga q se mueve en un campo magnético no uniforme

${\displaystyle \overrightarrow{B}=-\alpha yB\widehat{j}+(1+\alpha z)B\widehat{k}}$

La ecuación del movimiento es

${\displaystyle \begin{array}{l}m\frac{d^2\overrightarrow{r}}{d{t}^2}=q\overrightarrow{v}\times \overrightarrow{B}\\ \frac{d^2\overrightarrow{r}}{d{t}^2}=\frac{qB}{m}\left|\begin{array}{ccc}\widehat{i}& \widehat{j}& \widehat{k}\\ \frac{dx}{dt}& \frac{dy}{dt}& \frac{dz}{dt}\\ 0& -\alpha y& (1+\alpha z)\end{array}\right|\end{array}}$

Que da lugar a un sistema de tres ecuaciones diferenciales

${\displaystyle \{\begin{array}{l}\frac{d^{2}x}{d{t}^2}=\omega \left((1+\alpha z)\frac{dy}{dt}+\alpha y\frac{dz}{dt}\right)\\ \frac{d^{2}y}{d{t}^2}=-\omega (1+\alpha z)\frac{dx}{dt}\\ \frac{d^{2}z}{d{t}^2}=-\omega \alpha y\frac{dx}{dt}\end{array}}$

con ω=qB/m. La primera ecuación diferencial se puede integrar

${\displaystyle \begin{array}{l}\frac{d^{2}x}{d{t}^2}=\omega \frac{dy}{dt}+\omega \alpha \frac{d(y·z)}{dt}\\ \frac{dx}{dt}=\omega (1+\alpha z)y+C\end{array}}$

Donde C es una constante de integración que determinaremos a partir de las condiciones iniciales en el instante t=0

Conocido dx/dt, las otras dos ecuaciones diferenciales se expresan

${\displaystyle \{\begin{array}{l}\frac{d^{2}y}{d{t}^2}=-{\omega }^2{\left(1+\alpha z\right)}^{2}y-\omega C(1+\alpha z)\\ \frac{d^{2}z}{d{t}^2}=-{\omega }^2(1+\alpha z)\alpha y^2-\omega \alpha Cy\end{array}}$

Tenemos, por tanto, un sistema formado por una ecuación diferencial de primer orden y de dos ecuaciones diferenciales de segundo orden, que vamos a expresar de forma adimensional. Multiplicamos la primera ecuación diferencial por α/ω, la segunda y tercera por α/ω², y expresamos las ecuaciones diferenciales en términos de X=α·x, Y=α·y, Z=α·z, τ=ω·t y D=αC/ω

${\displaystyle \{\begin{array}{l}\frac{dX}{d\tau }=(1+Z)Y+D\\ \frac{d^{2}Y}{d{\tau }^2}=-{\left(1+Z\right)}^{2}Y-D(1+Z)\\ \frac{d^{2}Z}{d{\tau }^2}=-(1+Z)Y^2-DY\end{array}}$

Resolvemos este sistema de ecuaciones diferenciales con las siguientes condiciones iniciales: En el instante τ=0, la posición inicial es X=0, Y=Y₀, X=0, la velocidad inicial es dX/dτ=U_x, dY/dτ=0, dZ/dτ=U_z.

La constante de integración D que aprece en la primera ecuación diferencial, se determina a partir de las condiciones iniciales, en el instante τ=0, U_x=Y₀+D

Resultados

Sea Y₀=0.66, U_x=0.02, U_z=0.08. Resolvemos el sistema de ecuaciones diferenciales hasta el tiempo τ=140

• Coordenada X en función del tiempo τ

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot(t,x(:,1))
grid on
xlabel('\tau')
ylabel('X');
title('X en función de \tau')



Comprobamos que la energía de la partícula se mantiene constante, ya que la fuerza que ejerce el campo magnético $\overrightarrow{F}=q\overrightarrow{v}\times \overrightarrow{B}$ es siempre perpendicular al desplazamiento $\overrightarrow{dr}$ de la partícula

${\displaystyle {\left(\frac{dX}{d\tau }\right)}^2+{\left(\frac{dY}{d\tau }\right)}^2+{\left(\frac{dZ}{d\tau }\right)}^2=U_x^2+U_z^2}$

>> ((1+x(:,4)).*x(:,2)+D).^2+x(:,3).^2+x(:,5).^2
ans =
    0.0068
    0.0068      
    ......   
    0.0067
    0.0067
    0.0067
>> 0.02^2+0.08^2
ans =    0.0068 

• Coordenada Y en función del tiempo τ

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot(t,x(:,2))
grid on
xlabel('\tau')
ylabel('Y');
title('Y en función de \tau')



• Coordenada Z en función del tiempo τ

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot(t,x(:,4))
grid on
xlabel('\tau')
ylabel('Z');
title('Z en función de \tau')



La coordenada Z alcanza un máximo Z_m=2.33 en el tiempo τ_m=69.3

>> [xmax, nmax]=max(x(:,4))
xmax =    2.3292
nmax =   858
>> t(nmax)
ans =   69.2907

• La trayectoria en el espacio X, Y, Z

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot3(x(:,1),x(:,2),x(:,4))
grid on
xlabel('X')
ylabel('Y');
zlabel('Z')
title('Trayectoria X,Y,Z')
view(40,20)



• La trayectoria en el plano Z, Y

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot(x(:,4),x(:,2))
%plot(x(:,1),x(:,2))
grid on
xlabel('Z')
ylabel('Y');
title('Trayectoria Z,Y')



La partícula cargada se refleja en el instante τ_m en la posición Z_m

• La trayectoria en el plano X, Y

%X=x(1), Y=x(2), dY/dt=x(3), Z=x(4), dZ/dt=x(5)
x0=[0, 0.66, 0, 0, 0.08]; %condiciones iniciales
tf=140; %tiempo final
D=0.02-x0(2); %constante de integración
fg=@(t,x)[(1+x(4))*x(2)+D; x(3); -(1+x(4))^2*x(2)-D*(1+x(4)); 
x(5); -(1+x(4))*x(2)^2-D*x(2)];
[t,x]=ode45(fg,[0,tf],x0);
plot(x(:,1),x(:,2))
grid on
xlabel('X')
ylabel('Y');
title('Trayectoria X,Y')

Referencias

O. L. de Lange, J. Pierrus Solved Problems in Classical Mechanics. Analytical and numerical solutions with comments. Question 7.23, pp. 201-205