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

Eric Anterrieu

412 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

LabVIEW," *Int. J. Eng. Educ.*, vol 21, no 1, pp. 49-54.

http://www.mathworks.com/matlabcentral/fileexchange

vol 48, no 1, pp. 37–46.

pp. 173-184.

25.

[9] J. L. Guzmán, K. J. Åström, S. Dormido, T. Hägglund, and Y. Piguet, (2008), Interactive learning modules for PID control, *IEEE Cont. Syst. Mag.*, vol. 28, no. 5, pp. 118-134. [10] T. Kikuchi, T. Kenjo and S. Fukuda, (2002), "Developing an educational simulation program for the PM stepping motor," *IEEE Trans. Educ.*, vol 45, no 1, pp-70-78. [11] S. Ayasun and C. O. Nwankpa, (2005), Induction motor tests using Matlab/Simulink and their integration into undergraduate electric machinery courses, *IEEE Trans. Educ.*,

[12] N. Patrascoiu, (2005), Modeling and Simulation of the DC Motor Using Matlab and

[13] D. Hercog and K. Jezernik, (2005), Rapid Control Prototyping using Matlab/Simulink

[15] Casini, M. y D. Prattichizzo, (2003), The automatic control Telelab: A user-friendly

[16] R. Pastor, C. Martín, J. Sánchez, and S. Dormido, (2005), Development of a XML-based lab for remote control experiments on a servo motor, *Int. J. Elect. Eng. Educ.*, vol. 42, nº 2,

[17] R. Kelly and J. Moreno, (2001), Learning PID structures in an introductory course of

[21] J. Apkarian and K. J. Astrom, (2004), A laptop servo for control education, *IEEE Cont.* 

[24] B. M. Olds, B. M. Moskal and R. L. Miller, (2005), Assessment in Engineering Education: Evolution, Approaches, and Future Collaborations," *J. Eng. Educ.*, vol. 94, no. 1, pp. 13-

and DSP-based Motor Control," *Int. J. Eng. Educ*., vol. 21 no. 4, pp.596-605.

interface for distance learning," *IEEE Trans. Educ.*, vol.46, no. 2, pp. 252-257.

[18] LM628/LM629 Precision Motion Controller, [on-line] Retrieved March 2012 from:

[14] MATLAB-Central-File-Exchange [on-line], Retrieved March 2012 from:

automatic control," *IEEE Trans. Educ.*, vol. 44, no. 4, pp. 373-376.

[19] MAGELLAN Motion Control ICs, [on-line] Retrieved March 2012 from: http://www.pmdcorp.com/downloads/Magellan\_PUG\_v23\_Feb\_2008.pdf

[22] K. J. Åström and T. Hägglund, (2005), Advanced PID Control, ISA-Society. [23] The Matlab interactive module for servo systems learning: [on-line]:

http://cache.national.com/ds/LM/LM628.pdf

http://www.galilmc.com/about/index.html

*Syst. Mag.*, vol. 24, no. 5, pp-70-73.

[20] GALIL Controllers, [on-line] Retrieved March 2012 from:

http://www.esp.uem.es/aliane/servosystems/ss.zip

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46457

## **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*) as well as fully implicit ones of the form *f*(*t*, *y*, *y*� ) = 0. Generally

©2012 Anterrieu, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Anterrieu, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 415

The iconal equation [4], which describes the path of the light propagating in a medium with

where **r** denotes the position vector of a point on the ray and d*s* is an element of length along the path. From the point of view of the numerical implementation of a ray tracing procedure,

When light is propagating in a plane, the vector equation (2) reduces to a system of two scalar

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

with

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

d*x* <sup>d</sup> <sup>=</sup> *px*,

d*y* <sup>d</sup> <sup>=</sup> *py*.

= ∇*n*, (1)

<sup>∇</sup>*n*2. (2)

(3)

(4)

d d*s* � *n* d**r** d*s* �

> d2**r** <sup>d</sup><sup>2</sup> <sup>=</sup> <sup>1</sup> 2

equations of the second order. We therefore have in Cartesian coordinates:

⎧ ⎪⎪⎪⎪⎨ d2*x* <sup>d</sup><sup>2</sup> <sup>=</sup> <sup>1</sup> 2 *∂n*<sup>2</sup> *∂x* ,

d2*y* <sup>d</sup><sup>2</sup> <sup>=</sup> <sup>1</sup> 2 *∂n*<sup>2</sup> *∂y* .

⎪⎪⎪⎪⎩

cost.

**3.1. Solving the iconal equation**

refractive index *n*, is an ODE of the second order:

ODEs of the first order, (3) has to be rewritten:

⎧ ⎪⎪⎪⎨ d*px* <sup>d</sup> <sup>=</sup> <sup>1</sup> 2 *∂n*<sup>2</sup> *∂x* ,

d*py* <sup>d</sup> <sup>=</sup> <sup>1</sup> 2 *∂n*<sup>2</sup> *∂y* ,

⎪⎪⎪⎩

values of *px*, *py*, d*px*/d and d*py*/d for the same value of .

it is more convenient to set d = d*s*/*n* so that (1) now reads:

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:


The output arguments contain the solution approximated at discrete points:


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**

2 Will-be-set-by-IN-TECH

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

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,

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

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."

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,

equation has to be rewritten as an equivalent system of *n* first-order ODEs.

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

The output arguments contain the solution approximated at discrete points:

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

common to all the solver functions, is:

The input arguments are:

of the solver.

• t, a vector of time points.

**3. Non-linear optics**

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

• 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.

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{\mathbf{d}}{\mathbf{ds}} \left( n \frac{\mathbf{dr}}{\mathbf{ds}} \right) = \nabla n\_\prime \tag{1}$$

where **r** denotes the position vector of a point on the ray and d*s* 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 = d*s*/*n* so that (1) now reads:

$$\frac{\mathbf{d}^2 \mathbf{r}}{\mathbf{d}\ell^2} = \frac{1}{2} \nabla n^2. \tag{2}$$

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{cases} \frac{\text{d}^2 \text{x}}{\text{d}\ell^2} = \frac{1}{2} \frac{\partial n^2}{\partial \text{x}},\\\\ \frac{\text{d}^2 y}{\text{d}\ell^2} = \frac{1}{2} \frac{\partial n^2}{\partial y}. \end{cases} \tag{3}$$

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{cases} \frac{\mathrm{d}p\_{\mathrm{x}}}{\mathrm{d}\ell} = \frac{1}{2} \frac{\partial n^{2}}{\partial \mathbf{x}}, & \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\ell} = p\_{\mathrm{x}\prime} \\\\ \frac{\mathrm{d}p\_{\mathrm{y}}}{\mathrm{d}\ell} = \frac{1}{2} \frac{\partial n^{2}}{\partial \mathbf{y}}, & \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}\ell} = p\_{\mathrm{y}}. \end{cases} \tag{4}$$

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*, d*px*/d and d*py*/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*<sup>2</sup> 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*, (d*x*/d)*<sup>i</sup>* and (d*y*/d)*<sup>i</sup>* at a given starting point *Mi*. The values of the solution, namely *x*, *y*, d*x*/d et d*y*/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:

$$m^2(\rho) = 1 + \frac{\rho\_o^2}{\rho^2} \quad \text{with} \quad \rho = \sqrt{x^2 + y^2}. \tag{5}$$

so that no discontinuity occurs at the boundary. According to (5), within D the components

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup> <sup>2</sup>*xρ*<sup>2</sup>

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup> <sup>2</sup>*yρ*<sup>2</sup>

where Ro is used for storing the value of parameter *ρ<sup>o</sup>* and R that of the radius *R* of D.

and d*y* = d*s* sin *α*. Since at any point on the ray path d = d*s*/*n*, we therefore have:

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 d*x* = d*s* cos *α*

(d*x*/d)*<sup>i</sup>* = *ni* cos *α* and (d*y*/d)*<sup>i</sup>* = *ni* sin *α*.

Finally, the radius of D is set, for example, to *R* = 4*ρ<sup>o</sup>* with *ρ<sup>o</sup>* = 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*ρ<sup>o</sup>* from the *Ox* axis and with *xi* = 6*ρo*. The following lines are the MATLAB code corresponding to all these initial settings.

whereas they are null outside D since *n* is there uniform and equal to *ni* given by (6). The

0 (*x*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*2)<sup>2</sup> ,

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 417

0 (*x*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*2)<sup>2</sup> , (7)

*∂n*<sup>2</sup>

*∂n*<sup>2</sup>

of the gradient of *n*<sup>2</sup> in Cartesian coordinates are:

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

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

function dn2 is therefore reduced to:

if (sqrt(x^2+y^2) > R) dn2dx = 0; dn2dy = 0;

Ro = param(1); R = param(2);

else

end

>> 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];

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 D with radius *R*. Outside D, *n* is supposed to be uniform and equal to

*ni* = 1 + *ρ*<sup>2</sup> *<sup>o</sup>*/*R*2, (6)

so that no discontinuity occurs at the boundary. According to (5), within D the components of the gradient of *n*<sup>2</sup> in Cartesian coordinates are:

$$\begin{split} \frac{\partial n^2}{\partial x} &= -\frac{2x\rho\_0^2}{(\varkappa^2 + y^2)^2}, \\ \frac{\partial n^2}{\partial y} &= -\frac{2y\rho\_0^2}{(\varkappa^2 + y^2)^2}, \end{split} \tag{7}$$

whereas they are null outside D since *n* is there uniform and equal to *ni* 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
```
4 Will-be-set-by-IN-TECH

The students have to supply a function dn2 which should return the Cartesian components of the gradient of *n*<sup>2</sup> 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 solver is fed with a vector fi which contains the initial settings of the problem, namely *xi*, *yi*, (d*x*/d)*<sup>i</sup>* and (d*y*/d)*<sup>i</sup>* at a given starting point *Mi*. The values of the solution, namely *x*, *y*, d*x*/d et d*y*/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*.

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

*o*

*<sup>ρ</sup>*<sup>2</sup> with *<sup>ρ</sup>* <sup>=</sup>

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

Coming back to (5), the region where the media is non-homogeneous is restricted to a disk D

*x*<sup>2</sup> + *y*2. (5)

*<sup>o</sup>*/*R*2, (6)

*<sup>n</sup>*2(*ρ*) = <sup>1</sup> <sup>+</sup> *<sup>ρ</sup>*<sup>2</sup>

supplied by the students: this is exactly what is done in a classroom.

with radius *R*. Outside D, *n* is supposed to be uniform and equal to

*ni* = 1 + *ρ*<sup>2</sup>

function df=iconalODE(l,f) x = f(1); px = f(3); y = f(2); py = f(4); % gradient of n^2

dpxdl = 0.5\*dn2dx; dpydl = 0.5\*dn2dy;

>> X = f(:,1); >> Y = f(:,2);

% ode

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

the built-in function ode45 provided by MATLAB:

**3.2. Application in non-homogeneous media**

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

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

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

where Ro is used for storing the value of parameter *ρ<sup>o</sup>* and R that of the radius *R* of D.

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 d*x* = d*s* cos *α* and d*y* = d*s* sin *α*. Since at any point on the ray path d = d*s*/*n*, we therefore have:

(d*x*/d)*<sup>i</sup>* = *ni* cos *α* and (d*y*/d)*<sup>i</sup>* = *ni* sin *α*.

Finally, the radius of D is set, for example, to *R* = 4*ρ<sup>o</sup>* with *ρ<sup>o</sup>* = 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*ρ<sup>o</sup>* 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:

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

before leaving the non-homogeneous media.

0

situation is represented in Figure 4 for *h* = 0.965 *ρo*.

COULOMB potential [5].

/2

π

π

3 /2

π

Γ

2

π

5 /2

π

3

π

deviation is about 25◦ with respect to the incident direction.

**Figure 2**

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 D 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

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 419

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,

1 1.5 2 2.5 3 3.5 4

*h*/*ρo*

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 D. 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 D. Such an amazing

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

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

```
>> figure(1);
>> plot(X/Ro,Y/Ro,'r-');
```
**Figure 1.** An example of a ray path (in red) between points *Mi* and *Me*. Outside the disk D (in grey) the media is homogeneous with a constant refractive index: the light is travelling in straight line. On the contrary, the media in D is non-homogeneous with spatially varying refractive index: the path of the light is curved.

As expected, outside the disk D, the light is travelling in straight line since the media is homogeneous with uniform refractive index given by (6). On the contrary, within D 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 D, 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 D 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:

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

that is to say in MATLAB language:

418 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>7</sup> Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 419

```
>> gamma = atan(dYdX(end))+pi
```
6 Will-be-set-by-IN-TECH

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

−6 −4 −2 <sup>0</sup> <sup>2</sup> 4 6 −6

*x*/*ρo*

**Figure 1.** An example of a ray path (in red) between points *Mi* and *Me*. Outside the disk D (in grey) the media is homogeneous with a constant refractive index: the light is travelling in straight line. On the contrary, the media in D is non-homogeneous with spatially varying refractive index: the path of the

As expected, outside the disk D, the light is travelling in straight line since the media is homogeneous with uniform refractive index given by (6). On the contrary, within D 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 D, 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

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

Γ = arctan(d*y*/d*x*)*<sup>e</sup>* + *π*

with respect to the location where the distance from the center of D is minimal.

<sup>Γ</sup> ❄ *<sup>h</sup>*

*Mi* •

✻

Figure 1. It corresponds to the following lines:

**Figure 1**

−4

−2

0

*y*/*ρo*

*Me* •

2

4

6

>> figure(1);

light is curved.

media according to:

that is to say in MATLAB language:

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

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 D 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 D. 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 D. 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 4**

*h*/*ρ<sup>o</sup>* = 0.965

−4

**4. Relativistic electrodynamics**

−2

0

*y*/*ρo*

2

4

6

−6 −4 −2 <sup>0</sup> <sup>2</sup> 4 6 −6

*x*/*ρo*

**Figure 4.** Ray path for *h* = 0.965 *ρ<sup>o</sup>* (central part is enlarged on the top-left).

**Figure 4**

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 421

−0.6

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.

−0.4

−0.2

*y*/*ρo*

0

0.2

0.4

0.6

−0.6 −0.4 −0.2 0 0.2 0.4 0.6

*x*/*ρo*

8 Will-be-set-by-IN-TECH 420 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>9</sup>

**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 *ρ<sup>o</sup>* (central part is enlarged on the top-left).

## **4. Relativistic electrodynamics**

8 Will-be-set-by-IN-TECH

*h*/*ρ<sup>o</sup>* = 1.1195

*h*/*ρ<sup>o</sup>* = 1.0582 Γ = 3*π*/2

*h*/*ρ<sup>o</sup>* = 1.0019 Γ = 3*π*

Γ = *π*

**Figure 3.** Some light trajectories for various values of *h*. According to the inverse return of light, all the

**Figure 3**

*h*/*ρ<sup>o</sup>* = 1.0123 Γ = 5*π*/2

−6 −4 −2 <sup>0</sup> <sup>2</sup> 4 6 −6

*x*/*ρo*

trajectories are valid solutions for propagation in both directions.

−4

−2

0

*y*/*ρo*

2

4

6

*h*/*ρ<sup>o</sup>* = 1.2986 Γ = *π*/2

*h*/*ρ<sup>o</sup>* = 1.0289 Γ = 2*π*

> 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{d\mathbf{p}}{dt} = q(\mathbf{E} + \mathbf{v} \times \mathbf{B}),\tag{8}$$

function df=odeEB(t,f,q,m) % particle coordinates

P = sqrt(Px.^2+Py.^2+Pz.^2); E = sqrt(P^2\*c^2+m^2\*c^4);

% electric and magnetic fields [Ex,Ey,Ez]=fieldE(x,y,z,t); [Bx,By,Bz]=fieldB(x,y,z,t);

dPxdt = q\*(Ex + Vy\*Bz - Vz\*By); dPydt = q\*(Ey + Vz\*Bx - Vx\*Bz); dPzdt = q\*(Ez + Vx\*By - Vy\*Bx);

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

>> x = f(:,1); Px = f(:,4); >> y = f(:,2); Py = f(:,5); >> z = f(:,3); Pz = f(:,6);

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

odeEB is solved again with the aid of the built-in function ode45:

the momentum **p** at the corresponding time, are returned in the array f.

This function is general enough for solving various problems without re-writing any code since the mass *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

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 423

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*)*<sup>i</sup>* as well as the initial components of the momentum (*px*, *py*, *pz*)*<sup>i</sup>* 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

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*.

% particle momentum

% particle energy

% particle velocity Vx = Px\*c^2/E; Vy = Py\*c^2/E; Vz = Pz\*c^2/E;

x = f(1); y = f(2); z = f(3);

Px = f(4); Py = f(5); Pz = f(6);

% ode

where **v** is the velocity of the particle, **p** = *γ m***v** is the relativistic expression of the momentum with *<sup>γ</sup>* <sup>=</sup> 1/√<sup>1</sup> <sup>−</sup> **<sup>v</sup>**2/*c*<sup>2</sup> the relativistic factor and *<sup>c</sup>* 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{cases} \frac{\mathbf{d}p\_x}{\mathbf{d}t} = q(E\_x + v\_y B\_z - v\_z B\_y) \\\\ \frac{\mathbf{d}p\_y}{\mathbf{d}t} = q(E\_y + v\_z B\_x - v\_x B\_z) \\\\ \frac{\mathbf{d}p\_z}{\mathbf{d}t} = q(E\_z + v\_x B\_y - v\_y B\_x) \end{cases} \tag{9}$$

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 E of the free particle:

$$\mathbf{v} = \frac{\mathbf{p}}{\gamma m} = \frac{\mathbf{p}c^2}{\mathcal{E}} = \mathbf{p} \frac{c^2}{\sqrt{\mathbf{p}^2 c^2 + m^2 c^4}}.\tag{10}$$

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

$$\begin{cases} 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{cases} \tag{11}$$

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 d*x*/d*t*, d*y*/d*t* and d*z*/d*t* of the speed vector **v** as well as the components of the first derivative of the momentum d*px*/d*t*, d*py*/d*t* and d*pz*/d*t* at the same time.

422 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>11</sup> Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 423

```
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];
```
10 Will-be-set-by-IN-TECH

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

where **v** is the velocity of the particle, **p** = *γ m***v** is the relativistic expression of the momentum with *<sup>γ</sup>* <sup>=</sup> 1/√<sup>1</sup> <sup>−</sup> **<sup>v</sup>**2/*c*<sup>2</sup> the relativistic factor and *<sup>c</sup>* is the speed of light. In Cartesian coordinates, the previous vector equation can be put in the form of a system of three scalar

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> *<sup>q</sup>*(*Ex* <sup>+</sup> *vyBz* <sup>−</sup> *vzBy*)

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> *<sup>q</sup>*(*Ey* <sup>+</sup> *vzBx* <sup>−</sup> *vxBz*)

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> *<sup>q</sup>*(*Ez* <sup>+</sup> *vxBy* <sup>−</sup> *vyBx*)

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

<sup>E</sup> <sup>=</sup> **<sup>p</sup>**

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

� (*p*<sup>2</sup> *<sup>x</sup>* + *p*<sup>2</sup>

� (*p*<sup>2</sup> *<sup>x</sup>* + *p*<sup>2</sup>

� (*p*<sup>2</sup> *<sup>x</sup>* + *p*<sup>2</sup>

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 d*x*/d*t*, d*y*/d*t* and d*z*/d*t* of the speed vector **v** as well as the components of the first derivative of the

*c*2

*<sup>z</sup>* )*c*<sup>2</sup> + *m*2*c*<sup>4</sup>

*<sup>z</sup>* )*c*<sup>2</sup> + *m*2*c*<sup>4</sup>

*<sup>z</sup>* )*c*<sup>2</sup> + *m*2*c*<sup>4</sup>

*<sup>y</sup>* + *p*<sup>2</sup>

*<sup>y</sup>* + *p*<sup>2</sup>

*<sup>y</sup>* + *p*<sup>2</sup>

�**p**2*c*<sup>2</sup> <sup>+</sup> *<sup>m</sup>*2*c*<sup>4</sup> . (10)

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> *<sup>q</sup>*(**<sup>E</sup>** <sup>+</sup> **<sup>v</sup>** <sup>×</sup> **<sup>B</sup>**), (8)

(9)

(11)

d**p**

⎧ ⎪⎪⎪⎪⎪⎪⎪⎨ d*px*

d*py*

d*pz*

LORENTZ factor *γ*. This is done by introducing the energy E of the free particle:

*<sup>γ</sup> <sup>m</sup>* <sup>=</sup> **<sup>p</sup>***c*<sup>2</sup>

⎪⎪⎪⎪⎪⎪⎪⎩

**<sup>v</sup>** <sup>=</sup> **<sup>p</sup>**

*vx* = *pxc*2/

*vy* = *pyc*2/

*vz* = *pzc*2/

⎧ ⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎩

momentum d*px*/d*t*, d*py*/d*t* and d*pz*/d*t* at the same time.

**4.1. Solving the movement equation**

ODEs of the first order:

differential equation (ODE) of the first order [7]:

This function is general enough for solving various problems without re-writing any code since the mass *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*)*<sup>i</sup>* as well as the initial components of the momentum (*px*, *py*, *pz*)*<sup>i</sup>* 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.

>> Vi = 2\*c/3;

>> Pyi = 0; >> Pzi = 0; >> Xi = 0;

>> Zi = 0;

>> figure(5);

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

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

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

returned by the ode45 solver is represented on Figure 5:

**Figure 5**

*Mi*

−0.04

**E** = *E* **e***z*, **B** = *B* **e***<sup>z</sup>* and **v***<sup>i</sup>* = *vi***e***<sup>x</sup>* with *vi* = 2 *c*/3.

*z*

−10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0

−0.02

0 0.02

*x*

0.04 −0.04

**Figure 5.** Trajectory of a relativistic electron in an electromagnetic field when **E** and **B** are collinear:

to *Oy*) along the *Oz* axis (in the opposite direction of the **E** because here *q* < 0).

The trajectory is that of a circular helix with a variable step which rolls anticlockwise (from *Ox*

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

−0.02 0 0.02

*y*

0.04

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

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 425

## **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***<sup>z</sup>* and **B** = *B* **e***z*, and uniform with intensities *<sup>E</sup>* and *<sup>B</sup>* set to 3 · 106 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 *Mi* with an initial speed vector **v***<sup>i</sup>* = *vi***e***<sup>x</sup>* 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* = *ρ<sup>i</sup>* and *zi* = 0, for example, where *ρ<sup>i</sup>* = *γimvi*/*qB* is the initial radius of curvature. We therefore have the following initial conditions:

424 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>13</sup> Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 425

```
>> Vi = 2*c/3;
>> Pxi = gamma(Vi)*m*Vi;
>> Pyi = 0;
>> Pzi = 0;
>> Xi = 0;
>> Yi = Pxi/(q*0.05);
>> Zi = 0;
```
12 Will-be-set-by-IN-TECH

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,

The applications presented in this section are concerned with the movement of a relativistic

Depending on the orientation of **E** with respect to that of **B**, the trajectory of the electron may

In this first case, **E** and **B** are collinear, oriented along the *Oz* axis, **E** = *E* **e***<sup>z</sup>* and **B** = *B* **e***z*, and uniform with intensities *<sup>E</sup>* and *<sup>B</sup>* set to 3 · 106 V/m and 0.05 T, respectively. The functions

The electron is injected in the region where the electromagnetic field is applied from a point *Mi* with an initial speed vector **v***<sup>i</sup>* = *vi***e***<sup>x</sup>* 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* = *ρ<sup>i</sup>* and *zi* = 0, for example, where *ρ<sup>i</sup>* = *γimvi*/*qB* is the initial radius of curvature.

fieldE and fieldB supplied by the students are therefore reduced to:

electron under the simultaneous action of an electric field **E** and a magnetic one **B**.

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

**4.2. Applications in uniform fields**

>> c = 2.99792458E+08; % m/s >> q = -1.60217649E-19; % C >> m = 9.10938215E-31; % kg

function [Ex,Ey,Ez]=fieldE(x,y,z,t) % cartesian components of electric field

function [Bx,By,Bz]=fieldB(x,y,z,t) % cartesian components of magnetic field

We therefore have the following initial conditions:

>> Vx = Px\*c^2./E; >> Vy = Py\*c^2./E; >> Vz = Pz\*c^2./E;

Vy and Vz.

be very different.

Ex = 0; Ey = 0;

Bx = 0; By = 0; Bz = 0.05;

and

Ez = 3.00E+06;

*4.2.1.* **E** *and* **B** *collinear*

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:

```
>> 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***<sup>z</sup>* and **v***<sup>i</sup>* = *vi***e***<sup>x</sup>* with *vi* = 2 *c*/3.

The trajectory is that of a circular helix with a variable step which rolls anticlockwise (from *Ox* 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*<sup>2</sup> *<sup>x</sup>* + *v*<sup>2</sup> *<sup>y</sup>* + *v*<sup>2</sup> *<sup>z</sup>* )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*:

Bx = 0; By = 0.02; Bz = 0;

direction, namely −**e***<sup>x</sup>* because here *q* < 0.

−0.06

0

**E** = *E* **e***x*, **B** = *B* **e***<sup>y</sup>* and **v***<sup>i</sup>* = *vi***e***<sup>z</sup>* with here *vi* = 2 *c*/3.

0.5

1

1.5

*z*

2

2.5

−0.03

0 0.03

*x*

0.06 −0.06

*Mi*

**Figure 7.** Trajectory of a relativistic electron in an electromagnetic field when **E** and **B** are perpendicular:

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

−0.03 0 0.03

*y*

0.06

**Figure 7**

The electron is injected in the region where the electromagnetic field is applied from a point *Mi* with a initial speed vector **v***<sup>i</sup>* = *vi***e***<sup>z</sup>* 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 **v***i*, **E** and **B**. Indeed, only the magnetic force can periodically return the particle since *q* **v** × **B** operates in the 0*xz* 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

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 427

```
>> 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***<sup>x</sup>* and **B** = *B* **e***y*, and uniform with intensities *<sup>E</sup>* and *<sup>B</sup>* set to 3 · 106 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
```
426 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>15</sup> Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 427

> Bx = 0; By = 0.02; Bz = 0;

14 Will-be-set-by-IN-TECH

work on the particle. Consequently, as expected from the theory of relativity, the students can

trajectory increases rapidly (under the only action of the electric field **E**) from *vi* = 2 *c*/3 up to

0 5 10 15 20 25 30 35

*t* (ns)

Now **E** and **B** are perpendicular to each other, **E** = *E* **e***<sup>x</sup>* and **B** = *B* **e***y*, and uniform with intensities *<sup>E</sup>* and *<sup>B</sup>* set to 3 · 106 V/m and 0.02 T, respectively. The functions fieldE and

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

*<sup>x</sup>* + *v*<sup>2</sup>

*<sup>y</sup>* + *v*<sup>2</sup>

*<sup>z</sup>* )1/2 along the

easily observed on Figure 6 that the velocity of the particle *v* = (*v*<sup>2</sup>

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

**Figure 6**

0.7

fieldB are modified accordingly by the students:

function [Ex,Ey,Ez]=fieldE(x,y,z,t) % cartesian components of electric field

function [Bx,By,Bz]=fieldB(x,y,z,t) % cartesian components of magnetic field

*4.2.2.* **E** *and* **B** *perpendicular*

Ex = 3.00E+06;

Ey = 0; Ez = 0;

and

0.8

0.9

*v*/*c*

1

the very limit *c*:

>> figure(6);

The electron is injected in the region where the electromagnetic field is applied from a point *Mi* with a initial speed vector **v***<sup>i</sup>* = *vi***e***<sup>z</sup>* 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 **v***i*, **E** and **B**. Indeed, only the magnetic force can periodically return the particle since *q* **v** × **B** operates in the 0*xz* 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***<sup>x</sup>* 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***<sup>y</sup>* and **v***<sup>i</sup>* = *vi***e***<sup>z</sup>* with here *vi* = 2 *c*/3.

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 [12] which is a device used for filtering charged particles with a given velocity (or energy), for example in electron microscopes or in

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 429

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

Here, the particle is under the lonely action of a static but non-uniform magnetic field **B**, the

*xz z*2 *o*

*yz z*2 *o*

> *z*2 *z*2 *o* )

(12)

*<sup>x</sup>* + *v*<sup>2</sup>

*<sup>y</sup>* + *v*<sup>2</sup>

*<sup>x</sup>* + *v*<sup>2</sup>

*<sup>z</sup>* )1/2 along

*<sup>y</sup>*)1/2 of

*Bx* = −*Bo*

*By* = −*Bo*

*Bz* = *Bo*(1 +

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 **v***<sup>i</sup>* = *vi*(sin *αi***e***<sup>x</sup>* + cos *αi***e***z*) with a velocity *vi* = 2 *c*/3 and *α<sup>i</sup>* = 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

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

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

easily introduced. This is exactly what is done by the students in the classroom.

⎧ ⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎩

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 *v* = (*v*<sup>2</sup>

the trajectory as well as those of the axial and radial components <sup>|</sup>*vz*<sup>|</sup> and *<sup>v</sup><sup>ρ</sup>* = (*v*<sup>2</sup>

spectrometers.

components of which are:

the speed vector **v**:

not work.

>> 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-');

**4.3. Applications in time/space varying fields**

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

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**

16 Will-be-set-by-IN-TECH

**Figure 8**

−0.06

0

0.5

1

*z*

1.5

2

2.5

−0.03

*x*

0 0.03

identifies to that of a straight line.

0.06 −0.03

*vi* = 0.58 *c*

*vi* = 0.46 *c*

0

*y*

0.03

**Figure 8.** Same as Figure 7 for various initial velocities *vi*: when *vi* = *E*/*B*, the trajectory of the electron

*vi* = 0.54 *c*

*vi* = *<sup>E</sup>*/*B* 0.5 *c*

> 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{cases} B\_X = -B\_o \frac{\chi z}{z\_o^2} \\\\ B\_Y = -B\_o \frac{yz}{z\_o^2} \\\\ B\_z = B\_o(1 + \frac{z^2}{z\_o^2}) \end{cases} \tag{12}$$

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 **v***<sup>i</sup>* = *vi*(sin *αi***e***<sup>x</sup>* + cos *αi***e***z*) with a velocity *vi* = 2 *c*/3 and *α<sup>i</sup>* = 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 *v* = (*v*<sup>2</sup> *<sup>x</sup>* + *v*<sup>2</sup> *<sup>y</sup>* + *v*<sup>2</sup> *<sup>z</sup>* )1/2 along the trajectory as well as those of the axial and radial components <sup>|</sup>*vz*<sup>|</sup> and *<sup>v</sup><sup>ρ</sup>* = (*v*<sup>2</sup> *<sup>x</sup>* + *v*<sup>2</sup> *<sup>y</sup>*)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 *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

18 Will-be-set-by-IN-TECH 430 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Illustrating Amazing Effects of Modern Physics with Numerical Simulations Conducted in the Classroom <sup>19</sup>

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<sup>ρ</sup>* 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 [1]. It is encountered, for example, in a magnetic bottle which is the superposition of two magnetic

Illustrating Amazing Eff ects of Modern Physics with Numerical Simulations Conducted in the Classroom 431

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

[1] Allen R.J. (1962). A demonstration of the magnetic mirror effect, *American Journal of*

[2] Anterrieu E. & Pérez J.-P. (2011). How to make attractive the teaching of relativistic electrodynamics?, *Proceedings of IEEE International Conference on Computer as a Tool*,

[3] Anterrieu E. & Pérez J.-P. (2010). Illustrating amazing effects of optics with the computer, *Proceedings of IEEE Engineering Education Conference*, pages 180-185, Madrid (Spain),

[4] Born M. & Wolf E. (1999). *Principles of optics*, 7th edition, Cambridge University Press,

[5] Boyer T.H. (2004). Unfamiliar trajectories for a relativistic particle in a Kepler or Coulomb potential, *American Journal of Physics*, Vol. 72, No. 8, pages 992-997, 2004.

mirrors used together for confinement purpose.

and it may vary from one student to another.

*National de la Recherche Scientifique & Université de Toulouse, France*

*Physics*, Vol. 30, No. 12, pages 867-869, 1962.

pages 1-4, Lisbon (Portugal), 27-29 April 2011.

**5. Conclusion**

**Author details**

Eric Anterrieu

**6. References**

14-16 April 2010.

Cambridge.

**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.

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<sup>ρ</sup>* 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 [1]. It is encountered, for example, in a magnetic bottle which is the superposition of two magnetic mirrors used together for confinement purpose.

## **5. Conclusion**

18 Will-be-set-by-IN-TECH

**Figure 9**

−0.02

0.1

0.2

0.3

0.4

*v*/*c*

0.5

0.6

0.7

−3

−2

−1

0

*z*

1

2

3

−0.01

0 0.01

*x*

**Figure 10**

0.02 −0.02

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

<sup>0</sup> <sup>25</sup> <sup>50</sup> <sup>75</sup> <sup>100</sup> <sup>125</sup> <sup>150</sup> <sup>0</sup>

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

*t* (ns)

−0.01 0 0.01

*y*

0.02

*v vρ vz* 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.

## **Author details**

Eric Anterrieu *National de la Recherche Scientifique & Université de Toulouse, France*

## **6. References**

	- [6] Dormand J.R. & Prince P.J. (1980). A family of embedded Runge-Kutta formulae, *Journal of Computational and Applied Mathematics*, Vol. 6, No. 1, pages 19-26, 1980.
	- [7] Feynman R.P. (2005). *The Feynman lectures on physics: the definitive and extended edition*, 2nd edition, Addison Wesley, Boston.
	- [8] Hoefer W.J.R. (1973). On teaching relativity to undergraduate electrical engineers, *IEEE Transactions on Education*, Vol. 16, No. 3, pages 152-156, 1973.
	- [9] Luneburg R.K. (1964). *Mathematical theory of optics*, 1st edition, University of California Press, Berkeley.
	- [10] Petouris A., Humphries J., Russell P., Vourdas A. & Jones G.R (1994). Computer aided design for teaching modern optics to electronic engineering undergraduates, *IEE Proceedings Science, Measurement & Technology*, Vol. 141, No. 3, pages 171-176, 1994.
	- [11] Sader U. & Jodl H.J. (1987). Computer-aided physics-particles in electromagnetic fields, *European Journal of Physics*, Vol. 8, No. 2, pages 88-92, 1987.
	- [12] Wien W. (1898). Untersuchungen über die electrische entladung in verdünnten gasen, *Annalen der Physik*, Vol. 301, No. 6, pages 440-452, 1898.

© 2012 Lessa et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Lessa et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**The MatLabTM Software Application** 

**and Power System** 

http://dx.doi.org/10.5772/48555

**1. Introduction** 

**for Electrical Engineering Simulations** 

Leonardo da S. Lessa, Afonso J. Prado, Rodrigo Cleber Silva, Sérgio Kurokawa, Luiz F. Bovolato and José Pissolato Filho

Transmission lines are the elements of the electric power system that connect the load to the generation joining the production facilities of energy over large geographic areas. You could say that the transmission of electricity is one of the most important contributions that

The distribution of current potential differences and the transfer of energy along a transmission line can be analyzed by various methods, and it is expected that all lead to the same result. In engineering problems, in general, it can not indiscriminately apply a single formula for solving a specific problem, without the full knowledge of the limitations and simplifications accepted in its derivation. It is worth mentioning that this circumstance would lead to its misuse. The so-called mathematical solutions of physical phenomena

Therefore, there are several models that represent the transmission lines and can be classified as to the nature of their parameters in the model parameters constant and variable

The parameters in the models in terms of frequency are easy to use, but can not adequately represent the line in the entire range of frequencies at which these phenomena are transient in nature. In most cases, these designs increase the amplitude of the high order harmonics,

For adequate representation of the transmission line should be considered that the longitudinal line parameters are strongly frequency dependent, including the variable

Additional information is available at the end of the chapter

engineering has offered to the modern civilization.

typically require simplifications and idealizations.

distorting the waveform peaks and producing exaggerated errors.

parameters the models with frequency.
