Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom

1. Introduction

Real experiments are sometimes missing or highly expensive to implement when teaching modern physics. This is an important issue, especially when the theories are not intuitive. However, with the incorporation of computers into the classroom it has become much easier to illustrate some amazing effects of modern physics that cannot be brought to the attention of students without the aid of numerical simulations. Thanks to the performance of MATLAB langage, as well as to the easiness to code a theoretical problem, the students can really conduct numerical expriments and play with them almost in real time. This computer-based teaching approach is illustrate here with the capabilities to use MATLAB for solving initial values problems of ordinary differential equations encountered in relativistic electrodynamics [2] as well as in non-linear optics [3].

2. Ordinary differential equations in Matlab

In mathematics, an ordinary differential equation (or ODE) is a relation f that contains one or more derivatives of a dependent variable y with respect to a single independent variable t, usually referred to as time. Many studies have been devoted to the solution of ODEs. In the case where the equation is linear, it can be solved by analytical methods. Unfortunately, most of the interesting ODEs are non-linear and, with a few exceptions, cannot be solved analytically. However, in science and in engineering, a numeric approximation to the solution is often good enough to solve a problem. Numerical integration is that part of numerical analysis which studies the numerical solutions of ODEs with the aid of iterative solvers.

MATLAB offers several numerical algorithms to solve a wide variety of ODEs. These solvers are designed to handle only explicit ODEs of the form y’ = f (t, y), linearly implicit ODEs of the form M(t, y) y’ = f (t, y) aswell as fully implicit ones of the form f (t, y, y’) = 0. Generally there are many functions that satisfy a given ODE, and additional information is necessary to specify the solution of interest. In an initial value problem, the solution of interest satisfies a specific initial condition y(t_o) = y_o at a given initial time t_o. Finally, the ODE solvers available in MATLAB accept only first-order ODEs. However, ODEs often involve derivatives of order higher than one. As a consequence, to use the ODE solvers for solving an ODE of order n, the equation has to be rewritten as an equivalent system of n first-order ODEs.

All the ODE solvers proposed in MATLAB share a common syntax that makes them easy to try if it is not apparent which one is the most appropriate. To apply a different method to the same problem, simply change the name of the ODE solver function. The simplest syntax, common to all the solver functions, is:

>> [t,y] = solver(odefun,tspan,y0,options);

where solver is one of the ODE solver functions available in MATLAB.

The input arguments are:

• odefun, a handle to a function that evaluates the ODE.

• tspan, a vector specifying the interval of integration.

• y0, a vector of initial conditions for the problem.

• options, a structure of optional parameters to change the default integration properties of the solver.

The output arguments contain the solution approximated at discrete points:

• t, a vector of time points.

• y, a solution array with the values of the solution y at the time t returned in t.

Beginning with initial condition, the solver steps through the time interval, computing a solution at each time step. If the solution for a time step satisfies the solver’s error tolerance criteria, it is a successful step. Otherwise, it is a failed attempt; the solver shrinks the step size and tries again. The solver ode45 is based on a RUNGE-KUTTA formula [6] and is according to MATLAB documentation “the best function to apply as a first try for most problems.”

3. Non-linear optics

Ray optics is the branch of optics [4] in which all the subtle wave effects are neglected: the light is considered as travelling along rays which can only change their direction by refraction or reflection. The usefullness of the computer for studying, or designing, optical systems within an educational frame, is well known and it has been already suggested [10], but only in the paraxial approximation of geometrical optics. Indeed, when light propagates in media with constant refractive index, the SNELL-DESCARTES laws can be applied for implementing a fast numerical ray-tracing procedure based on a geometrical approach of the problem. However, when propagation takes place in media with non-homogeneous refractive index, the differential equation governing the propagation of light has to be solved with the aid of the computer. This section describes a fast, accurate and easy to use code for illustrating the propagation of light in such media, sometimes with amazing effects that cannot be brought to the attention of students without the aid of numerical simulations, or except at an expensive cost.

3.1. Solving the iconal equation

The iconal equation [4], which describes the path of the light propagating in a medium with refractive index n, is an ODE of the second order:

$\frac{\text{d}}{\text{d}s}\left(n\frac{\text{d}r}{\text{d}s}\right)=\nabla n,$E1

where r denotes the position vector of a point on the ray and ds is an element of length along the path. From the point of view of the numerical implementation of a ray tracing procedure, it is more convenient to set dℓ = ds/n so that (1) now reads:

$\frac{{\text{d}}^{2}r}{\text{d}{\ell }^2}=\frac{1}{2}\nabla n^2.$E2

When light is propagating in a plane, the vector equation (2) reduces to a system of two scalar equations of the second order. We therefore have in Cartesian coordinates:

$\{\begin{array}{ccc}\frac{{\text{d}}^{2}x}{\text{d}{\ell }^2}& =& \frac{1}{2}\frac{\partial n^2}{\partial x},\\ \frac{{\text{d}}^{2}y}{\text{d}{\ell }^2}& =& \frac{1}{2}\frac{\partial n^2}{\partial y}.\end{array}$E3

As soon as the refractive index does not vary with x and y, (3) is straiforward to integrate and leads to a straight line path. When n is not uniform and according to the dependency of n with respect to x and y, it may not be possible to integrate (3) without the aid of the computer. This is why the capabilities of MATLAB to quickly solve an ODE is used here for studying the propagation of light in non-homogeneous media. However, since MATLAB can only solve ODEs of the first order, (3) has to be rewritten:

$\{\begin{array}{ccc}\frac{\text{d}{p}_x}{\text{d}\ell }& =& \frac{1}{2}\frac{\partial n^2}{\partial x},\\ \frac{\text{d}{p}_y}{\text{d}\ell }& =& \frac{1}{2}\frac{\partial n^2}{\partial y},\end{array}\text{with}\begin{array}{ccc}\frac{\text{d}x}{\text{d}\ell }& =& p_x,\\ \frac{\text{d}y}{\text{d}\ell }& =& p_y.\end{array}$E4

The MATLAB function iconalODE contains the final definition of the problem: in vector f are stored the values of x, y, p_x and p_y for a given value of ℓ, whereas the vector df returns the values of p_x, p_y, dp_x/d ℓ and dp_y/d ℓ for the same value of ℓ.

function df=iconalODE(l,f)

x = f(1); px = f(3);

y = f(2); py = f(4);

% gradient of n^2

[dn2dx, dn2dy] = dn2(x,y,param);

% ode

dpxdl = 0.5*dn2dx;

dpydl = 0.5*dn2dy;

df = [px; py; dpxdl; dpydl];

The students have to supply a function dn2 which should return the Cartesian components of the gradient of n² at any point (x, y), the expression of which may depend on parameters stored in the param vector. The equation described in iconalODE is solved with the aid of the built-in function ode45 provided by MATLAB:

>> [l,f] = ode45(’iconalODE’,l,fi);

>> X = f(:,1);

>> Y = f(:,2);

>> dYdX = f(:,4)./f(:,3);

The solver is fed with a vector fi which contains the initial settings of the problem, namely x_i, y_i, (dx/dℓ)_i and (dy/dℓ)_i at a given starting point M_i. The values of the solution, namely x, y, dx/dℓ et dy/dℓ along the ray path, are returned in the array f: here, for convenience reasons, vectors X, Y and dYdX are used for storing the values of x, y and dy/dx.

3.2. Application in non-homogeneous media

In this study, the refractive index of the non-homogeneous media satisfies the law:

$n^2(\rho )=1+\frac{{\rho }_o^2}{{\rho }^2}\text{with}\rho =\sqrt{x^2+y^2}.$E5

Another choice would have been to study the propagation of light in LUNEBURG lens [9], however contrary to expression (5) the refractive index would have not present a singularity and consequently the example would have been less constraining from the computational point of view and less illustrative from the point of view of the propagation of light in non-homogeneous media. However, according to the approach adopted here, moving from one study to another just requires to change the lines of code written in the function dn2 supplied by the students: this is exactly what is done in a classroom.

Coming back to (5), the region where the media is non-homogeneous is restricted to a disk with radius R. Outside , n is supposed to be uniform and equal to

$n_i=\sqrt{1+{\rho }_o^2/R^2},$E6

so that no discontinuity occurs at the boundary. According to (5), within the components of the gradient of n² in Cartesian coordinates are:

$\begin{array}{ccc}\frac{\partial n^2}{\partial x}& =& -\frac{2x{\rho }_0^2}{{(x^2+y^2)}^2},\\ \frac{\partial n^2}{\partial y}& =& -\frac{2y{\rho }_0^2}{{(x^2+y^2)}^2},\end{array}$E7

whereas they are null outside since n is there uniform and equal to n_i given by (6). The function dn2 is therefore reduced to:

function [dn2dx,dn2dy]=dn2(x,y,param)

Ro = param(1);

R = param(2);

if (sqrt(x^2+y^2) > R)

dn2dx = 0;

dn2dy = 0;

else

dn2dx = -2*x*Ro^2/(x^2+y^2)^2;

dn2dy = -2*y*Ro^2/(x^2+y^2)^2;

end

where Ro is used for storing the value of parameter ρ_o and R that of the radius R of .

Before solving the ODE described in functions iconalODE and dn2 it is necessary to set the initial conditions of the problem. Let us consider an incident ray starting at a point M_i(x_i, y_i) and making an angle α with the Ox axis. In the neighbourhood of M_i, we have dx = ds cos α and dy = ds sin α. Since at any point on the ray path dℓ = ds/n, we therefore have:

${\left(\text{d}x/\text{d}\ell \right)}_i=n_i\cos \alpha \text{and}{\left(\text{d}y/\text{d}\ell \right)}_i=n_i\sin \alpha .$E8

Finally, the radius of is set, for example, to R = 4ρ_o with ρ_o = 1. Since the system exhibits a central symmetry, the study is restricted here to incidents rays parallel to the Ox axis, so that α = π, and starting from M_i located at a distance y_i = h = 2ρ_o from the Ox axis and with x_i = 6ρ_o. The following lines are the MATLAB code corresponding to all these initial settings.

>> Ro = 1;

>> R = 4*Ro;

>> h = 2*Ro;

>> Xi = 6*Ro;

>> Yi = h;

>> alpha = pi;

>> Ni = sqrt(1 + (Ro/R)^2);

>> dXidl = Ni*cos(alpha);

>> dYidl = Ni*sin(alpha);

>> fi = [Xi Yi dXidl dYidl];

After point M_i, the light travels according to the ODE described in the functions iconalODE and dn2. The complete path of the ray returned by the ode45 solver is represented on Figure 1. It corresponds to the following lines:

>> figure(1);

>> plot(X/Ro,Y/Ro,’r-’);

Figure 1.

An example of a ray path (in red) between points M_i and M_e. Outside the disk (in grey) the media is homogeneous with a constant refractive index: the light is travelling in straight line. On the contrary, the media in is non-homogeneous with spatially varying refractive index: the path of the light is curved.

As expected, outside the disk , the light is travelling in straight line since the media is homogeneous with uniform refractive index given by (6). On the contrary, within the path is curved since the light is here travelling in a non-homogeneous media with spatially varying refractive index according to (5). The curvature of the path is oriented towards the center of the disk , that is to say towards the direction of the gradient of n. The ray continuously tends towards the center without reaching it. According to the inverse return of light, the path of the light obtained when solving the iconal equation (1) from M_i to M_e and that obtained when solving it from M_e to M_i in the reverse direction of propagation are nothing but the same. In addition, owing to the central symmetry of the system, the path of the light is symmetrical with respect to the location where the distance from the center of is minimal.

Coming back to the situation of Figure 1, the deviation angle Γ after travelling from M_i through the non-homogeneous media can be computed at any point M_e in the homogeneous media according to:

$$\Gamma =\arctan {(\text{d}y/\text{d}x)}_e+\pi$$E9

that is to say in MATLAB language:

>> gamma = atan(dYdX(end))+pi

As observed by the students, h is playing the role of an impact parameter with respect to the singularity of the refractive index at the center of the disk and the angle Γ that of a scattering angle with regards to the direction of the incoming ray. Shown in Figure 1, for h = 2ρ_o, the deviation is about 25° with respect to the incident direction.

An interesting and convenient aspect of conducting numerical simulations with the aid of the computer in a classroom with MATLAB is the capability to re-run the code with different initial settings and to obtain the new path of the light almost immediately. This is exactly what is done by the students for investigating the influence of the parameter h on the deviation angle Γ. As shown in Figure 2, Γ varies with h in a non-linear manner. Values of Γ greater than π/2 correspond to rays which make a half turn, or even more than a complete turn, before leaving the non-homogeneous media.

Figure 2.

Variations of the deviation (or scattering) angle Γ with the (impact) parameter h.

Some particular situations are shown in Figure 3 for different values of the ratio h/ρ_o. Although trajectories for large values of Γ are not without surprise for the students, they usually accept that the light can make a large number of turns before leaving the disk . However, their reaction is simply incredible when they discover, even with numerical experiences, that below a limit for the ratio h/ρ_o the path of the light identifies to that of a spiral: the ray seems to be attracted by the singularity of the refractive index of the non-homgeneous media for ρ = 0 and does not emerge from the disk . Such an amazing situation is represented in Figure 4 for h = 0.965 ρ_o.

Coming back to the analogy between h and an impact parameter and between Γ and a scattering angle, these trajectories remind those of a particule captured by a KEPLER or COULOMB potential [5].

Figure 3.

Some light trajectories for various values of h. According to the inverse return of light, all the trajectories are valid solutions for propagation in both directions.

Figure 4.

Ray path for h = 0.965 ρ_o (central part is enlarged on the top-left).

4. Relativistic electrodynamics

Classical electrodynamics is an important branch of the theoretical physics [7], which has numerous applications in physical engineering. On one hand, the usefulness of the computer for studying the movement of particles in electromagnetic fields is now well known and it has been already suggested [11], but only in the Newtonian approximation. However, the trajectory becomes complicated to establish when the electric and magnetic fields are acting together, whatever the initial velocity. On the other hand, the teaching of relativity to undergraduate electrical engineers appears also to be a necessity [8]. A stimulated way to gain the attention and the interest of students for relativity, which is not an intuitive theory, is precisely to begin by studying the movement of fast particles in electromagnetic fields because the applications are numerous and surprising. This section describes a fast, accurate and easy to use code for computing the trajectories of charged relativistic particles under the action of the LORENTZ force and illustrating some effects of relativity that cannot be brought to the attention of students without the aid of numerical simulations, or except at an expensive cost.

4.1. Solving the movement equation

The equation which describes the trajectory of a relativistic particle of mass m and charge q under the simultaneous actions of an electric field E and a magnetic one B is an ordinary differential equation (ODE) of the first order [7]:

$\frac{\text{d}p}{\text{d}t}=q(E+v\times B),$E10

where v is the velocity of the particle, p = γmv is the relativistic expression of the momentum with $\gamma =1/\sqrt{1-v^2/c^2}$ the relativistic factor and c is the speed of light. In Cartesian coordinates, the previous vector equation can be put in the form of a system of three scalar ODEs of the first order:

$$\{\begin{array}{ccc}\frac{\text{d}{p}_x}{\text{d}t}& =& q(E_x+v_y{B}_z-v_z{B}_y)\\ \frac{\text{d}{p}_y}{\text{d}t}& =& q(E_y+v_z{B}_x-v_x{B}_z)\\ \frac{\text{d}{p}_z}{\text{d}t}& =& q(E_z+v_x{B}_y-v_y{B}_x)\end{array}$$E11

Owing to the factor γ, solving this system of ODEs is harder in the relativistic case than in the Newtonian one. Indeed, expression of the first derivatives of the components p_x, p_y et p_z of the momentum p is not simple because of the expression of the first derivative of γ itself. It is necessary to express the speed v as a function of the momentum p without involving the LORENTZ factor γ. This is done by introducing the energy ε of the free particle:

$v=\frac{p}{\gamma m}=\frac{p{c}^2}{ℰ}=p\frac{c^2}{\sqrt{p^2{c}^2+m^2{c}^4}}.$E12

This vector equation reduces again to three scalar equations in Cartesian coordinates:

$\{\begin{array}{ccc}{v}_x& =& p_x{c}^2/\sqrt{(p_x^2+p_y^2+p_z^2)c^2+m^2{c}^4}\\ v_y& =& p_y{c}^2/\sqrt{(p_x^2+p_y^2+p_z^2)c^2+m^2{c}^4}\\ v_z& =& p_z{c}^2/\sqrt{(p_x^2+p_y^2+p_z^2)c^2+m^2{c}^4}\end{array}$E13

Finally, solving (8) in Cartesian coordinates in the relativistic case is equivalent to solve the system of ODEs (9) together with (11). However, owing to the coupling between these equations, it may not be possible to integrate them without the aid of the computer. This is why the capabilities of MATLAB to quickly solve ODEs is used here for tackling this problem.

The MATLAB function odeEB contains the definition of the system to solve: in vector f are stored the coordinates (x, y, z) of the particle as well as the components (p_x, p_y, p_z) of the momentum p at a given time t, whereas the vector df returns the components dx/dt, dy/dt and dz/dt of the speed vector v as well as the components of the first derivative of the momentum dp_x/dt, dp_y/dt and dp_z/dt at the same time.

function df=odeEB(t,f,q,m)

% particle coordinates

x = f(1);

y = f(2);

z = f(3);

% particle momentum

Px = f(4);

Py = f(5);

Pz = f(6);

% particle energy

P = sqrt(Px.^2+Py.^2+Pz.^2);

E = sqrt(P^2*c^2+m^2*c^4);

% particle velocity

Vx = Px*c^2/E;

Vy = Py*c^2/E;

Vz = Pz*c^2/E;

% electric and magnetic fields

[Ex,Ey,Ez]=fieldE(x,y,z,t);

[Bx,By,Bz]=fieldB(x,y,z,t);

% ode

dPxdt = q*(Ex + Vy*Bz - Vz*By);

dPydt = q*(Ey + Vz*Bx - Vx*Bz);

dPzdt = q*(Ez + Vx*By - Vy*Bx);

df = [Vx; Vy; Vz; dPxdt; dPydt; dPzdt];

This function is general enough for solving various problems without re-writing any code since themass m and the electric charge q of the particle are given to the function. The students have only to supply the two functions fieldE and fieldB which should return the Cartesian components of E and B at any point (x, y, z) and at any time t. The equation described in odeEB is solved again with the aid of the built-in function ode45:

>> [t,f] = ode45(’odeEB’,t,fi);

The solver is fed with a vector fi which contains the initial settings of the problem, namely the initial position of the particle (x, y, z)_i as well as the initial components of the momentum (p_x, p_y, p_z)_i at a given starting point M_i. The values of the solution, namely the successive positions of the particle (x, y, z) along the trajectory as well as the components (p_x, p_y, p_z) of the momentum p at the corresponding time, are returned in the array f.

>> x = f(:,1); Px = f(:,4);

>> y = f(:,2); Py = f(:,5);

>> z = f(:,3); Pz = f(:,6);

Here, for convenience reasons, vectors x, y and z are used for storing the values of x, y and z. Likewise, vectors Px, Py and Pz are used for storing the values of p_x, p_y and p_z.

>> P = sqrt(Px.^2+Py.^2+Pz.^2);

>> E = sqrt(P.^2*c^2+m^2*c^4);

>> Vx = Px*c^2./E;

>> Vy = Py*c^2./E;

>> Vz = Pz*c^2./E;

The final step performed outside the solver is the computation of the components (v_x, v_y, v_z) of the speed vector v, according to (10). The corresponding values are written in vectors Vx, Vy and Vz.

4.2. Applications in uniform fields

The applications presented in this section are concerned with the movement of a relativistic electron under the simultaneous action of an electric field E and a magnetic one B.

>> c = 2.99792458E+08; % m/s

>> q = -1.60217649E-19; % C

>> m = 9.10938215E-31; % kg

Depending on the orientation of E with respect to that of B, the trajectory of the electron may be very different.

4.2.1. E and B collinear

In this first case, E and B are collinear, oriented along the Oz axis, E = E e_z and B = B e_z, and uniform with intensities E and B set to 3 10⁶ V/m and 0.05 T, respectively. The functions fieldE and fieldB supplied by the students are therefore reduced to:

function [Ex,Ey,Ez]=fieldE(x,y,z,t)

% cartesian components of electric field

Ex = 0;

Ey = 0;

Ez = 3.00E+06;

and

function [Bx,By,Bz]=fieldB(x,y,z,t)

% cartesian components of magnetic field

Bx = 0;

By = 0;

Bz = 0.05;

The electron is injected in the region where the electromagnetic field is applied from a point M_i with an initial speed vector v_i = v_ie_x perpendicular to E and B and a velocity v_i = 2 c/3. The coordinates of M_i are chosen in such a way that the trajectory rolls up along the Oz axis: x_i = 0, y_i = ρ_i and z_i = 0, for example, where ρ_i = γ_imv_i/qB is the initial radius of curvature. We therefore have the following initial conditions:

>> Vi = 2*c/3;

>> Pxi = gamma(Vi)*m*Vi;

>> Pyi = 0;

>> Pzi = 0;

>> Xi = 0;

>> Yi = Pxi/(q*0.05);

>> Zi = 0;

After point M_i, the electron is under the influence of E and B: its trajectory is the solution of the ODE described in function odeEB with the previous initial conditions. The complete path returned by the ode45 solver is represented on Figure 5:

>> figure(5);

>> plot3(x,y,z,’r-’);

Figure 5.

Trajectory of a relativistic electron in an electromagnetic field when E and B are collinear: E = E e_z, B = B e_z and v_i = v_ie_x with v_i = 2 c/3.

The trajectory is that of a circular helix with a variable step which rolls anticlockwise (fromOx to Oy) along the Oz axis (in the opposite direction of the E because here q < 0).

According to the LORENTZ force law (8), a charged particle can gain (or lose) energy from an electric field E, but not from a magnetic field B. Indeed, by definition the magnetic force q v × B is always perpendicular to the particle’s direction of motion and therefore does not work on the particle. Consequently, as expected from the theory of relativity, the students can easily observed on Figure 6 that the velocity of the particle $v={(v_x^2+v_y^2+v_z^2)}^{1/2}$ along the trajectory increases rapidly (under the only action of the electric field E) from vi = 2 c/3 up to the very limit c:

>> figure(6);

>> plot(t/1E-09,sqrt(Vx.^2+Vy.^2+Vz.^2)/c,’r-’);

Figure 6.

Relative velocity v/c of a relativistic electron along the trajectory shown in Figure 5.

4.2.2. E and B perpendicular

Now E and B are perpendicular to each other, E = E e_x and B = B e_y, and uniform with intensities E and B set to 3 10⁶ V/m and 0.02 T, respectively. The functions fieldE and fieldB are modified accordingly by the students:

function [Ex,Ey,Ez]=fieldE(x,y,z,t)

% cartesian components of electric field

Ex = 3.00E+06;

Ey = 0;

Ez = 0;

and

function [Bx,By,Bz]=fieldB(x,y,z,t)

% cartesian components of magnetic field

Bx = 0;

By = 0.02;

Bz = 0;

The electron is injected in the region where the electromagnetic field is applied from a point M_i with a initial speed vector v_i = v_ie_z perpendicular to both E and B and a velocity v_i = 2 c/3. The complete path returned by the ode45 solver is represented on Figure 7. The trajectory is now confined to a plan perpendicular to B and containing both E and the speed vector v. More exactly it is a sinusoidal curve with amplitude and period depending on v_i, E and B. Indeed, only the magnetic force can periodically return the particle since q v × B operates in the 0xz plan sometimes in a direction, sometimes in the other one, depending on the direction of v. On the other hand, the electric force q E also operates in this plan but always in the same direction, namely −e_x because here q < 0.

Figure 7.

Trajectory of a relativistic electron in an electromagnetic field when E and B are perpendicular: E = E e_x, B = B e_y and v_i = v_ie_z with here v_i = 2 c/3.

Figure 8.

Same as Figure 7 for various initial velocities v_i: when v_i = E/B, the trajectory of the electron identifies to that of a straight line.

An interesting and convenient aspect of conducting numerical simulations with the aid of the computer in a classroom with MATLAB is the capability to re-run the code with different initial settings and to obtain the new trajectory almost immediately. As observed by the students in such a situation, as soon as the value of v_i becomes closer to the ratio E/B, the amplitude of the oscillations decreases as illustrated in Figure 8. When v_i = E/B, the actions of the electric and magnetic forces offset each other so that the trajectory reduces to a straight line. Such an amazing situation is encountered in a WIEN filter [12] which is a device used for filtering charged particles with a given velocity (or energy), for example in electron microscopes or in spectrometers.

4.3. Applications in time/space varying fields

In the previous applications, the students dealt with uniform electromagnetic fields. Since the time and space coordinates are given to the functions fieldE and fieldB, some time and/or space variations in the intensity E and B of the electric and magnetic fields can be easily introduced. This is exactly what is done by the students in the classroom.

Here, the particle is under the lonely action of a static but non-uniform magnetic field B, the components of which are:

$\{\begin{array}{ccc}{B}_x& =& -B_o\frac{xz}{z_o^2}\\ B_y& =& -B_o\frac{yz}{z_o^2}\\ B_z& =& B_o(1+\frac{z^2}{z_o^2})\end{array}$E14

where B_o and z_o are taken equal to 0.05 T and 1 m, respectively. The functions fieldE and fieldB are modified accordingly by the students to reflect this new experimental setup.

The initial conditions of the problem are now the following ones. The electron is injected in the region where the magnetic field is applied from a point M_i with an initial speed vector v_i = v_i(sin α_ie_x + cos α_ie_z) with a velocity v_i = 2 c/3 and α_i = 20°, for example. The coordinates of M_i are chosen in such a way that the trajectory rolls along the Oz axis. The complete path returned by the ode45 solver is shown in Figure 9. The trajectory reminds that of a circular helix with a fixed step but with a variable radius of curvature ρ which raises its maximum for z = 0 where B_z raises its minimum and where B_x = B_y = 0.

Shown in Figure 10 are the variations of the velocity of the particle $v={(v_x^2+v_y^2+v_z^2)}^{1/2}$along the trajectory as well as those of the axial and radial components $\left|v_z\right|\text{and }{v}_{\rho }={(v_x^2+v_y^2)}^{1/2}$ of the speed vector v:

>> figure(10); hold(’on’);

>> plot(t/1E-09,abs(Vz)/c,’b-’);

>> plot(t/1E-09,sqrt(Vx.^2+Vy.^2)/c,’g-’);

>> plot(t/1E-09,sqrt(Vx.^2+Vy.^2+Vz.^2)/c,’r-’);

Contrary to what is observed in Figure 6 where the particule is moving under the simultaneous action of an electric field E and a magnetic one B, here the velocity v of the particle is constant and equal to v_i since it is under the lonely action of a (non-uniform) magnetic field B which cannot give (or take) energy to (or from) the particle because it does not work.

Coming back to the trajectory observed in Figure 9 in the light of Figure 10, as soon as B_z grows from z = 0, the particle draws circles around the Oz axis which become continually smaller and closer. At the same time, the axial contribution to the kinetic energy is converted into a rotational one up to v_z is cancelled out and v_ρ reach its maximum(v remaining constant and equal to v_i). Then, the particle turns and goes back, drawing again circles in the opposite direction which become larger and more distant up to z = 0. The particle has been reflected by the magnetic field B. Such an experimental setup is called a magnetic mirror [1]. It is encountered, for example, in a magnetic bottle which is the superposition of two magnetic mirrors used together for confinement purpose.

Figure 9.

Trajectory of a relativistic electron in the non-uniform magnetic field described in (12).

Figure 10.

Relative velocity v/c of a relativistic electron along the trajectory shown in Figure 9.

5. Conclusion

This contribution has described MATLAB codes for solving ordinary differential equations that cannot be solved without the aid of the computer in the field of non-linear optics as well as in relativistic electrodynamics. Since these two domains of modern physics are associated to non intuitive theories, concrete illustrations are welcomed to gain the attention and the interest of the students. The code is fast, accurate and easy to use. It has to be fed with few functions supplied by the students for switching from one study to another, the initial conditions of the problem being changed accordingly. Since the results are obtained almost in real time, these initial settings can be changed at will so that an interactive study is often conceivable with the computers in a classroom. Numerical examples have demonstrated the capabilities of these codes for illustrating in teaching conditions some amazing effects of modern physics that cannot be brought to the attention of students without the aid of the computer, or at an expensive cost. In a classroom setting, the time to be allocated between the underground physics, the modeling issues and the freewheel experimentation is of course left to the teacher and it may vary from one student to another.