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  as well as in non-linear optics .
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(to) = yo at a given initial time to. 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  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  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 , 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 , which describes the path of the light propagating in a medium with refractive index n, is an ODE of the second order:
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:
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:
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:
The MATLAB function iconalODE contains the final definition of the problem: in vector f are stored the values of x, y, px and py for a given value of ℓ, whereas the vector df returns the values of px, py, dpx/d ℓ and dpy/d ℓ for the same value of ℓ.
x = f(1); px = f(3);
y = f(2); py = f(4);
% gradient of n^2
[dn2dx, dn2dy] = dn2(x,y,param);
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 n2 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 xi, yi, (dx/dℓ)i and (dy/dℓ)i at a given starting point Mi. 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:
Another choice would have been to study the propagation of light in LUNEBURG lens , 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
so that no discontinuity occurs at the boundary. According to (5), within the components of the gradient of n2 in Cartesian coordinates are:
whereas they are null outside since n is there uniform and equal to ni given by (6). The function dn2 is therefore reduced to:
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 Mi(xi, yi) and making an angle α with the Ox axis. In the neighbourhood of Mi, we have dx = ds cos α and dy = ds sin α. Since at any point on the ray path dℓ = ds/n, we therefore have:
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 Mi located at a distance yi = h = 2ρo from the Ox axis and with xi = 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 Mi, 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:
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 Mi to Me and that obtained when solving it from Me to Mi 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 Mi through the non-homogeneous media can be computed at any point Me in the homogeneous media according to:
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.
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 .
4. Relativistic electrodynamics
Classical electrodynamics is an important branch of the theoretical physics , 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 , 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 . 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 :
where v is the velocity of the particle, p = γmv is the relativistic expression of the momentum with 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:
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 px, py et pz 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:
This vector equation reduces again to three scalar equations in Cartesian coordinates:
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 (px, py, pz) 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 dpx/dt, dpy/dt and dpz/dt at the same time.
% 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
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 (px, py, pz)i at a given starting point Mi. The values of the solution, namely the successive positions of the particle (x, y, z) along the trajectory as well as the components (px, py, pz) 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 px, py and pz.
>> 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 (vx, vy, vz) 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 ez and B = B ez, and uniform with intensities E and B set to 3 106 V/m and 0.05 T, respectively. The functions fieldE and fieldB supplied by the students are therefore reduced to:
% cartesian components of electric field
Ex = 0;
Ey = 0;
Ez = 3.00E+06;
% 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 Mi with an initial speed vector vi = viex perpendicular to E and B and a velocity vi = 2 c/3. The coordinates of Mi are chosen in such a way that the trajectory rolls up along the Oz axis: xi = 0, yi = ρi and zi = 0, for example, where ρi = γimvi/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 Mi, 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:
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 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:
4.2.2. E and B perpendicular
Now E and B are perpendicular to each other, E = E ex and B = B ey, and uniform with intensities E and B set to 3 106 V/m and 0.02 T, respectively. The functions fieldE and fieldB are modified accordingly by the students:
% cartesian components of electric field
Ex = 3.00E+06;
Ey = 0;
Ez = 0;
% 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 Mi with a initial speed vector vi = viez perpendicular to both E and B and a velocity vi = 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 vi, 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 −ex because here q < 0.
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 vi becomes closer to the ratio E/B, the amplitude of the oscillations decreases as illustrated in Figure 8. When vi = 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  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:
where Bo and zo 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 Mi with an initial speed vector vi = vi(sin αiex + cos αiez) with a velocity vi = 2 c/3 and αi = 20°, for example. The coordinates of Mi 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 Bz raises its minimum and where Bx = By = 0.
Shown in Figure 10 are the variations of the velocity of the particle along the trajectory as well as those of the axial and radial components of the speed vector v:
>> figure(10); hold(’on’);
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 vi 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 Bz 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 vz is cancelled out and vρ reach its maximum(v remaining constant and equal to vi). 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 . It is encountered, for example, in a magnetic bottle which is the superposition of two magnetic mirrors used together for confinement purpose.
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.