Accueil > Équipes > S2I > Eric Anterrieu > Ouvrages > Simulations numériques en relativité

 

Simulations numériques en relativité

Par Eric Anterrieu - 14/12/2005

 

AVIS IMPORTANT
Les programmes présentés ici sont destinés à illustrer les simulations numériques de l’annexe 3 du livre :
Relativité et invariance : fondements et applications
de José-Philippe PÉREZ avec la collaboration de Éric ANTERRIEU, DUNOD éditeur (réf. ISBN : 2-10-049173-3).

L’utilisation de ces programmes, écrits pour le logiciel de calcul scientifique MATLAB, est libre de droit et n’engage nullement la responsabilité des auteurs ni celle de l’éditeur ; elle ne peut en aucun cas se faire à des fins commerciales sans l’accord des auteurs et de l’éditeur, ce qui n’exclut donc pas d’utiliser ces programmes en classe ou en séance de travaux pratiques ou dirigés, à la condition expresse que la licence MATLAB soit à jour.

I.- Structures de données et fonctions de base

On construit d’abord des fonctions de base pour manipuler les facteurs relativistes :


function BETAe=beta(Ve)
% BETAe=BETA(Ve) returns the relativistic beta for a given velocity Ve.
%  » Ve = velocity in m/s
%  « BETAe = relativistic beta
%
% © Eric Anterrieu (2005)
%
% See also: INVBETA, GAMMA, INVGAMMA.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
BETAe = Ve/invbeta(1);
%
return;
beta.m


function Ve=invbeta(BETAe)
% Ve=INVBETA(BETAe) returns the velocity Ve for a given relativistic beta.
%  » BETAe = relativistic beta
%  « Ve = velocity in m/s
%
% © Eric Anterrieu (2005)
%
% See also: BETA, GAMMA, INVGAMMA.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
Ve = BETAe*299792458;
%
return;
invbeta.m


function GAMMAe=gamma(Ve)
% GAMMAe=GAMMA(Ve) returns the relativistic gamma for a given velocity Ve.
%  » Ve = velocity in m/s
%  « GAMMAe = relativistic gamma
%
% © Eric Anterrieu (2005)
%
% See also: INVGAMMA, BETA, INVBETA.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
BETAe = beta(Ve);
if (BETAe > 1-eps)
  GAMMAe = inf;
else
  GAMMAe = 1/sqrt(1-BETAe^2);
end
%
return;
gamma.m


function Ve=invgamma(GAMMAe)
% Ve=INVGAMMA(GAMMAe) returns the velocity Ve for a given relativistic gamma.
%  » GAMMAe = relativistic gamma
%  « Ve = velocity in m/s
%
% © Eric Anterrieu (2005)
%
% See also: GAMMA, BETA, INVBETA.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
if isinf(GAMMAe)
  Ve = invbeta(1);
else
  BETAe = sqrt(1-1/GAMMAe^2);
  Ve = invbeta(BETAe);
end
%
return;
invgamma.m

I.1.- Éléments constitutifs

On utilise les capacités de programmation orienté objet (POO) de MATLAB afin de mettre en oeuvre des classes destinées à manipuler les quadrivecteurs espace-temps et vitesse.

a) La classe QR

Le fichier QR.m dans le dossier @QR contient la définition de la classe QR pour manipuler des quadrivecteurs espace-temps :


function R4=QR(varargin)
% QR is the constructor function for QR object.
% The QR class is primarily meant to be used to store space-time quadrivector components.
% R4=QR(varargin)
% Called with one argument that is a QR object, the object is simply returned.
% Called with two arguments (ct,x), the constructor returns a new QR object.
% Called with three arguments (ct,x,y), the constructor returns a new QR object.
% Called with four arguments (ct,x,y,z), the constructor returns a new QR object.
%  » ct = time coordinate [m]
%  » x  = x coordinate [m]
%  » y  = y coordinate [m]
%  » z  = z coordinate [m]
%  « R4 = QR object
%
% © Eric Anterrieu (2004)
%
% See also: QR/DISPLAY.
%
c = invbeta(1);
%
switch (nargin)
case 1
  if (isa(varargin{1},'QR'))
    R4 = varargin{1};
  else
    error('Wrong argument type.');
  end
case 2
  R4.ct = c*varargin{1};
  R4.x  =   varargin{2};
  R4.y  = 0;
  R4.z  = 0;
  R4 = class(R4,'QR');
case 3
  R4.ct = c*varargin{1};
  R4.x  =   varargin{2};
  R4.y  =   varargin{3};
  R4.z  = 0;
  R4 = class(R4,'QR');
case 4
  R4.ct = c*varargin{1};
  R4.x  =   varargin{2};
  R4.y  =   varargin{3};
  R4.z  =   varargin{4};
  R4 = class(R4,'QR');
otherwise
  error('Wrong number of input arguments.');
 end
%
return;
@QR/QR.m

Dans le dossier @QR, on trouve la sur-définition de fonctions pour afficher/modifier les valeurs des propriétés d’un objet de la classe QR...


function display(R4)
% DISPLAY(R4) displays the components of a space-time quadrivector.
%  » R4 = quadrivector QR
%
% © Eric Anterrieu (2004)
%
% See also: QR/QR.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
disp(sprintf('ct = %g m',R4.ct));
disp(sprintf(' x = %g m',R4.x ));
disp(sprintf(' y = %g m',R4.y ));
disp(sprintf(' z = %g m',R4.z ));
%
return;
@QR/display.m


%
return;
@dioptre/set.m

Dans le dossier @QR, on trouve aussi des fonctions internes à la classe QR...


function s2=norm2(varargin)
% mat=MATRIX(obj) returns the refraction matrix of a dioptre.
%  » obj = dioptre object
%  « mat = refraction matrix
%
% © Eric Anterrieu (2004)
%
% See also: QR/QR, QR/GET.
%
switch (nargin)
case 1
if (isa(varargin{1},'QR'))
R4 = varargin{1};
else
error('Wrong argument type.');
end
s2 = (R4.ct)^2 - (R4.x)^2 - (R4.y)^2 - (R4.z)^2;
case 2
if (isa(varargin{1},'QR'))
R41 = varargin{1};
else
error('Wrong argument type.');
end
if (isa(varargin{2},'QR'))
R42 = varargin{2};
else
error('Wrong argument type.');
end
s2 = (R42.ct-R41.ct)^2 - (R42.x-R41.x)^2 - (R42.y-R41.y)^2 - (R42.z-R41.z)^2;
otherwise
error('Wrong number of input arguments.');
end
%
return;
@QR/norm2.m


function R42=lorentz(R41,Ve)
% R42=LORENTZ(R41,Ve) returns .
%  » R41 = 
%  » Ve = 
%  « R42 = 
%
% © Eric Anterrieu (2004)
%
% See also: QR/QR, QR/norm2.
%
BETAe =  beta(Ve);
GAMAe = gamma(Ve);
if isinf(GAMAe)
  error('Ve should be less than c');
end
%
R42.ct = GAMAe*(R41.ct + BETAe*R41.x );
R42.x  = GAMAe*(R41.x  + BETAe*R41.ct);
R42.y  = R41.y;
R42.z  = R41.z;
%
R42 = class(R42,'QR');
%
return;
@QR/lorentz.m

b) La classe QV

Le fichier QV.m dans le dossier @QV contient la définition de la classe QV pour manipuler des quadrivecteurs vitesse :


function V4=QV(varargin)
% QV is the constructor function for QV object.
% The QV class is primarily meant to be used to store spead quadrivector components.
% V4=QV(varargin)
% Called with one argument that is a QV object, the object is simply returned.
% Called with one argument (Vx), the constructor returns a new QV object.
% Called with two arguments (Vx,Vy), the constructor returns a new QV object.
% Called with three arguments (Vx,Vy,Vz), the constructor returns a new QV object.
%  » Vx = Vx component [m/s]
%  » Vy = Vy component [m/s]
%  » Vz = Vz component [m/s]
%  « V4 = QV object
%
% © Eric Anterrieu (2004)
%
% See also: QV/DISPLAY.
%
switch (nargin)
case 1
  if (isa(varargin{1},'QV'))
    V4 = varargin{1};
  else
    V4.gama = gamma(sqrt((varargin{1})^2));
    V4.Vx   = varargin{1};
    V4.Vy   = 0;
    V4.Vz   = 0;
    V4 = class(V4,'QV');
  end
case 2
  V4.gama = gamma(sqrt((varargin{1})^2+(varargin{2})^2));
  V4.Vx   = varargin{1};
  V4.Vy   = varargin{2};
  V4.Vz   = 0;
  V4 = class(V4,'QV');
case 3
  V4.gama = gamma(sqrt((varargin{1})^2+(varargin{2})^2+(varargin{3})^2));
  V4.Vx   = varargin{1};
  V4.Vy   = varargin{2};
  V4.Vz   = varargin{3};
  V4 = class(V4,'QV');
otherwise
  error('Wrong number of input arguments.');
 end
%
return;
@QV/QV.m

Dans le dossier @QV, on trouve la sur-définition de fonctions pour afficher/modifier les valeurs des propriétés d’un objet de la classe QV...


function display(V4)
% DISPLAY(V4) displays the components of a spead quadrivector.
%  » V4 = quadrivector QV
%
% © Eric Anterrieu (2004)
%
% See also: QV/QV.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
c = invbeta(1);
%
disp(sprintf('beta = %g',sqrt((V4.Vx)^2+(V4.Vy)^2+(V4.Vz)^2)/c));
disp(sprintf('Vx/c = %g',V4.Vx/c));
disp(sprintf('Vy/c = %g',V4.Vy/c));
disp(sprintf('Vz/c = %g',V4.Vz/c));
%
return;
@QV/display.m


function V42=lorentz(V41,Ve)
% V42=LORENTZ(V41,Ve) returns .
%  » V41 = 
%  » Ve = 
%  « V42 = 
%
% © Eric Anterrieu (2004)
%
% See also: QV/QV, QV/norm2.
%
c = invbeta(1);
BETAe =  beta(Ve);
GAMAe = gamma(Ve);
if isinf(GAMAe)
  error('Ve should be less than c');
end
%
GAMA1 = V41.gama;
if isinf(GAMA1)
  V42.gama = V41.gama;
else
  V42.gama = GAMAe*GAMA1*(1 + BETAe*V41.Vx/c);
end
V42.Vx = (Ve + V41.Vx)/(1 + BETAe*V41.Vx/c);
V42.Vy =      (V41.Vy)/(1 + BETAe*V41.Vx/c)/GAMAe;
V42.Vz =      (V41.Vz)/(1 + BETAe*V41.Vx/c)/GAMAe;
%
V42 = class(V42,'QV');
%
return;
@QV/lorentz.m

II.- Quadrivecteurs et transformation de Lorentz

A l’aide des éléments précédents, on illustre dans cette partie les propriétés des quadrivecteurs espace-temps et vitesse, relativement à la transformation de Lorentz. On vérifie aussi les relations liant la durée impropre à la durée propre Tp et le longueur impropre à la longueur propre Lp.

II.1.- Quadrivecteur espace-temps

Le script suivant calcul le carré de la pseudo norme séparant deux évènements A et B, faisant apparaitre d’abord un intervalle sAB de genre temps puis de genre espace, selon les coordonnées spatiaux-temporelles des évènements A et B :


c = invbeta(1);
%
A1 = QR(0.4/c,1)
B1 = QR(1.6/c,2)
SAB = norm2(A1,B1)
%
A1 = QR(0.4/c,1)
B1 = QR(1.2/c,2)
SAB = norm2(A1,B1)
%
return;
scriptII1.m

II.2.- Durée propre

Lorsque l’intervalle entre les deux évènements A et B est du genre temps, le script suivant calcule la durée propre Tp :


c = invbeta(1);
%
A1 = QR(0.4/c,1)
B1 = QR(1.6/c,2)
SAB = norm2(A1,B1)
%
cTp = Tpropre(A1,B1)
%
Ve = Rpropre(A1,B1);
A2 = lorentz(A1,-Ve)
B2 = lorentz(B1,-Ve)
%
return;
scriptII2.m

II.3.- Longueur propre

Le script suivant calcul la longueur propre Lp qui sépare deux évènements A et B :


c = invbeta(1);
%
A1 = QR(0.4/c,1)
B1 = QR(1.6/c,2)
SAB = norm2(A1,B1)
%
A1 = QR(0.4/c,1)
B1 = QR(1.2/c,2)
SAB = norm2(A1,B1)
%
Lp = Lpropre(A1,B1)
%
Ve = Rpropre(A1,B1);
A2 = lorentz(A1,-Ve)
B2 = lorentz(B1,-Ve)
%
return;
scriptII3.m

II.4.- Quadrivecteur vitesse

Le script suivant


c = invbeta(1);
%
Ue = invbeta(0.5);
U1 = invbeta(0.8);
V41 = QV(0,U1)
V42 = lorentz(V41,Ue)
%
PN1 = sqrt(norm2((V41))/c
PN2 = sqrt(norm2((V42))/c
clear V41
V41 = lorentz(V42,-Ue)
%
return;
scriptII4.m

III.- Emission de particules par un objet en mouvement

On s’intéresse à un noyau radioactif N, en mouvement rectiligne uniforme par rapport au référentiel du laboratoire, à la vitesse ue=ueex. Ce noyau émet, de manière isotrope, des particules P animées d’une vitesse v1 dans le référentiel galiléen R1 d’origine N. On se propose de calculer la distribution angulaire des vitesses v2 de P dans le référentiel R2 du laboratoire.

Le script suivant traite le cas relativiste où ue=0,5c et v1=0,75c, par exemple.


Ue = invbeta(0.5);
V1 = invbeta(0.75); % 0.25 / 0.5 / 0.75
T1 = (0:1:180)*pi/180;
for k=1:length(T1)
  V41=QV(V1*cos(T1(k)),V1*sin(T1(k)));
  V42=lorentz(V41,Ue);
  [T2(k),V2(k)]=polar(V42);
end
% variations of T2 with T1
figure(1); clf;
plot(T1*180/pi,T2*180/pi); axis([0 180 0 180]);
% variations of V2 with T1
figure(2); clf;
plot(T1*180/pi,V2); axis([0 180 0 1]);
%
return;
scriptIII.m

III.1.- Cas ultra-relativiste

Le script suivant traite le cas ultra-relativiste où ue=0,9c et v1=c.


Ue = invbeta(0.9);
V1 = invbeta(1);
T1 = (0:1:180)*pi/180;
for k=1:length(T1)
  V41=QV(V1*cos(T1(k)),V1*sin(T1(k)));
  V42=lorentz(V41,Ue);
  [T2(k),V2(k)]=polar(V42);
end
% variations of T2 with T1
figure(1); clf;
plot(T1*180/pi,T2*180/pi); axis([0 180 0 180]);
% variations of V2 with T1
figure(2); clf;
plot(T1*180/pi,V2); axis([0 180 0 1]);
%
return;
scriptIII1.m

III.2.- Approximation newtonienne

Le script suivant restitue l’approximation newtonienne avec ve=20 m/s et v1=30 m/s, par exemple.


Ve = 20;
V1 = 30; % 10 / 20 / 30
T1 = (0:1:180)*pi/180;
for k=1:length(T1)
  V41=QV(V1*cos(T1(k)),V1*sin(T1(k)));
  V42=lorentz(V41,Ve);
  [T2(k),V2(k)]=polar(V42);
end
c=invbeta(1);
% variations of T2 with T1
figure(1); clf;
plot(T1*180/pi,T2*180/pi); axis([0 180 0 180]);
% variations of V2 with T1
figure(2); clf;
plot(T1*180/pi,V2*c); axis([0 180 0 50]);
%
return;
scriptIII2.m

IV.- Particule dans un champ magnétique stationnaire

On s’intéresse ici au mouvement d’une particule de masse m et de charge q, soumise à la seule action d’un champ magnétique B stationnaire.

IV.1.- Particule dans un champ magnétique stationnaire uniforme

On se propose ici de résoudre l’équation différentielle vectorielle décrivant la trajectoire d’une particule, de masse m et de charge q, soumise à la seule action d’un champ magnétique B stationnaire uniforme. Le fichier odeB.m définit le problème que l’on cherche à résoudre, pour lequel on suppose que le champ magnétique B est orienté selon l’axe Oz.


function df=odeB(t,f,B,q,m)
% f = [x, y, z, Vx, Vy, Vz]
x = f(1); Vx = f(4);
y = f(2); Vy = f(5);
z = f(3); Vz = f(6);
% V and omega
V = sqrt(Vx^2 + Vy^2 + Vz^2);
w = (q*B)/(gamma(V)*m);
% EDO: dVx/dt, dVy/dt, dVz/dt
dVxdt =  w*Vy;
dVydt = -w*Vx;
dVzdt =  0;
% df = [Vx, Vy, Vz, dVx/dt, dVy/dt, dVz/dt]
df = [Vx; Vy; Vz; dVxdt; dVydt; dVzdt];
%
return;
odeB.m

Le script suivant précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier odeB.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire de la particule :


c  = invbeta(1);
% B
B = 0.05;
% q and m
q = -1.6021892E-19;
m =  0.9109534E-30;
% Vx, Vy and Vz at starting point
Vi = invbeta(2/3); alphai = 30*pi/180;
Vxi = Vi*sin(alphai);
Vyi = 0;
Vzi = Vi*cos(alphai);
% x, y and z at starting point
Zi = -2.5; Bi = B; Ri = (gamma(Vi)*m*Vi)/(q*Bi);
Yi = Ri;
Xi = 0;
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[t,f] = ode45('odeB',[0:0.01:9]/c,[Xi Yi Zi Vxi Vyi Vzi],options,B,q,m);
X = f(:,1); Vx = f(:,4);
Y = f(:,2); Vy = f(:,5);
Z = f(:,3); Vz = f(:,6);
% variations of [X, Y, Z] with ct
figure(1); clf;
subplot(3,1,1); plot(c*t,X,'b-');
subplot(3,1,2); plot(c*t,Y,'b-');
subplot(3,1,3); plot(c*t,Z,'b-');
% variations of [Vx, Vy, Vz]/c with ct
figure(2); clf;
subplot(3,1,1); plot(c*t,Vx/c,'b-');
subplot(3,1,2); plot(c*t,Vy/c,'b-');
subplot(3,1,3); plot(c*t,Vz/c,'b-');
% trajectory
figure(3); clf; hold('on');
plot3(X,Y,Z,'r-');
plot3([0 0],[0 0],[-2.5 2.5],'k--');
view(55,20); set(gca,'DataAspectRatio',[1 1 100]);
axis([-0.02 0.02 -0.02 0.02 -2.5 2.5]);
xlabel('X'); ylabel('Y'); zlabel('Z');
%
return;
scriptIV1.m

IV.2.- Focalisation dans un champ magnétique

.


c  = invbeta(1);
%
return;
scriptIV2.m

IV.3.- Focalisation dans un prisme magnétique

.


c  = invbeta(1);
%
return;
scriptIV3.m

IV.4.- Particule dans un champ magnétique stationnaire non uniforme

Comme dans le cas du champ stationnaire uniforme, on suppose que le champ magnétique B est orienté selon l’axe Oz : B=Bez. Mais ici, on suppose que son amplitude B varie spatialement selon la loi : B=Bo[1+(z/zo)2].
Le fichier odeB.m définit le problème que l’on cherche ainsi à résoudre :


function df=odeB(t,f,Bo,zo,q,m)
% f = [x, y, z, Vx, Vy, Vz]
x = f(1); Vx = f(4);
y = f(2); Vy = f(5);
z = f(3); Vz = f(6);
% B
B = Bo*(1+(z/zo)^2);
% V and omega
V = sqrt(Vx^2 + Vy^2 + Vz^2);
w = (q*B)/(gamma(V)*m);
% EDO: dVx/dt, dVy/dt, dVz/dt
dVxdt =  w*Vy;
dVydt = -w*Vx;
dVzdt =  0;
% df = [Vx, Vy, Vz, dVx/dt, dVy/dt, dVz/dt]
df = [Vx; Vy; Vz; dVxdt; dVydt; dVzdt];
%
return;
odeB.m

Le script suivant précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier odeB.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire de la particule :


c  = invbeta(1);
% B
Bo = 0.05;
zo = 1;
% q and m
q = -1.6021892E-19;
m =  0.9109534E-30;
% Vx, Vy and Vz at starting point
Vi = invbeta(2/3); alphai = 30*pi/180;
Vxi = Vi*sin(alphai);
Vyi = 0;
Vzi = Vi*cos(alphai);
% x, y and z at starting point
Zi = -2.5; Bi = Bo*(1+(Zi/zo).^2); Ri = (gamma(Vi)*m*Vi)/(q*Bi);
Yi = Ri;
Xi = 0;
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[t,f] = ode45('odeB',[0:0.01:9]/c,[Xi Yi Zi Vxi Vyi Vzi],options,Bo,zo,q,m);
X = f(:,1); Vx = f(:,4);
Y = f(:,2); Vy = f(:,5);
Z = f(:,3); Vz = f(:,6);
% variations of [X, Y, Z] with ct
figure(1); clf;
subplot(3,1,1); plot(c*t,X,'b-');
subplot(3,1,2); plot(c*t,Y,'b-');
subplot(3,1,3); plot(c*t,Z,'b-');
% variations of [Vx, Vy, Vz]/c with ct
figure(2); clf;
subplot(3,1,1); plot(c*t,Vx/c,'b-');
subplot(3,1,2); plot(c*t,Vy/c,'b-');
subplot(3,1,3); plot(c*t,Vz/c,'b-');
% trajectory
figure(3); clf; hold('on');
plot3(X,Y,Z,'r-');
plot3([0 0],[0 0],[-2.5 2.5],'k--');
view(55,20); set(gca,'DataAspectRatio',[1 1 100]);
axis([-0.02 0.02 -0.02 0.02 -2.5 2.5]);
xlabel('X'); ylabel('Y'); zlabel('Z');
%
return;
scriptIV4.m

V.- Action de champs électrique et magnétique

On se propose ici de résoudre l’équation différentielle vectorielle décrivant la trajectoire d’une particule, de masse m et de charge q, soumise à l’action simultanée d’un champ électrique E et d’un champ magnétique B.
Le fichier odeEB.m définit le problème que l’on cherche à résoudre, de la manière la plus générale puisqu’aucune hypothèse n’est faite sur l’orientation des champs E et B ni sur la nature de la particule :


function df=odeEB(t,f,E,B,q,m)
% f = [x, y, z, Px, Py, Pz]
x = f(1); Px = f(4);
y = f(2); Py = f(5);
z = f(3); Pz = f(6);
% E = [Ex, Ey, Ez] and B = [Bx, By, Bz]
Ex = E(1); Bx = B(1);
Ey = E(2); By = B(2);
Ez = E(3); Bz = B(3);
% V = [Vx, Vy, Vz]
c = invbeta(1);
Vx = Px*c^2/sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vy = Py*c^2/sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vz = Pz*c^2/sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
% EDO: dPx/dt, dPy/dt, dPz/dt
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, dPx/dt, dPy/dt, dPz/dt]
df = [Vx; Vy; Vz; dPxdt; dPydt; dPzdt];
%
return;
odeEB.m

V.1.- E et B colinéaires

Le script suivant étudie le cas où les champs électrique et magnétique sont colinéaires :
E=Eez et B=Bez avec E=0,03.108 V/m et B=0,05 T.
Il précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier odeEB.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire de la particule :


c  = invbeta(1);
% E and B
E = [0 0 1]*0.03E+08;
B = [0 0 1]*0.05E;
% q and m
q = -1.6021892E-19;
m =  0.9109534E-30;
% Px, Py and Pz at starting point
Vi = invbeta(2/3);
Pxi = gamma(Vi)*m*Vi;
Pyi = 0;
Pzi = 0;
% x, y and z at starting point
Ri = Pxi/(q*B(3));
Xi = 0;
Yi = Ri;
Zi = 0;
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[t,f] = ode45('odeEB',[0:0.01:10]/c,[Xi Yi Zi Pxi Pyi Pzi],options,E,B,q,m);
X = f(:,1); Px = f(:,4);
Y = f(:,2); Py = f(:,5);
Z = f(:,3); Pz = f(:,6);
Vx = Px*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vy = Py*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vz = Pz*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
% variations of [X, Y, Z] with ct
figure(1); clf;
subplot(3,1,1); plot(c*t,X,'b-');
subplot(3,1,2); plot(c*t,Y,'b-');
subplot(3,1,3); plot(c*t,Z,'b-');
% variations of [Vx, Vy, Vz]/c with ct
figure(2); clf;
subplot(3,1,1); plot(c*t,Vx/c,'b-');
subplot(3,1,2); plot(c*t,Vy/c,'b-');
subplot(3,1,3); plot(c*t,Vz/c,'b-');
% trajectory
figure(3); clf; hold('on');
plot3(X,Y,Z,'r-');
plot3([0 0],[0 0],[0 -10],'k--');
view(55,20); set(gca,'DataAspectRatio',[1 1 100]);
axis([-0.04 0.04 -0.04 0.04 -10 0]);
xlabel('X'); ylabel('Y'); zlabel('Z');
%
return;
scriptV1.m

V.2.- E et B perpendiculaires

Le script suivant étudie le cas où les champs électrique et magnétique sont perpendiculaires :
E=Eex et B=Bey avec E=0,03.108 V/m et B=0,02 T.
Il précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier odeEB.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire de la particule :


c  = invbeta(1);
% E and B
E = [1 0 0]*0.03E+08;
B = [0 1 0]*0.02E;
% q and m
q = -1.6021892E-19;
m =  0.9109534E-30;
% Px, Py and Pz at starting point
Vi = invbeta(2/3);
Pxi = 0;
Pyi = 0;
Pzi = gamma(Vi)*m*Vi;
% x, y and z at starting point
Xi = 0;
Yi = 0;
Zi = 0;
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[t,f] = ode45('odeEB',[0:0.01:5]/c,[Xi Yi Zi Pxi Pyi Pzi],options,E,B,q,m);
X = f(:,1); Px = f(:,4);
Y = f(:,2); Py = f(:,5);
Z = f(:,3); Pz = f(:,6);
Vx = Px*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vy = Py*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
Vz = Pz*c^2./sqrt((Px.^2+Py.^2+Pz.^2)*c^2+m^2*c^4);
% variations of [X, Y, Z] with ct
figure(1); clf;
subplot(3,1,1); plot(c*t,X,'b-');
subplot(3,1,2); plot(c*t,Y,'b-');
subplot(3,1,3); plot(c*t,Z,'b-');
% variations of [Vx, Vy, Vz]/c with ct
figure(2); clf;
subplot(3,1,1); plot(c*t,Vx/c,'b-');
subplot(3,1,2); plot(c*t,Vy/c,'b-');
subplot(3,1,3); plot(c*t,Vz/c,'b-');
% trajectory
figure(3); clf; hold('on');
plot3(X,Y,Z,'r-');
plot3([0 0],[0 0],[2.5 0],'k--');
view(55,20); set(gca,'DataAspectRatio',[1 1 100/6]);
axis([-0.06 0.06 -0.06 0.06 0 2.5]);
xlabel('X'); ylabel('Y'); zlabel('Z');
%
return;
scriptV2.m

 

 

[Page Précédente] [Dans la même rubrique] [Sommaire]

 

 

Site du CNRS Logo OMP
Logo UPS
Logo INSU-CNRS

Annuaire

Rechercher

Sur ce site



Bibliotheque

Laboratoire Astrophysique de Toulouse - Tarbes (UMR5572)

CNRS (Midi-Pyrénées)

Univ. Paul Sabatier

Liens utiles

Accueil Imprimer Contact mail Plan du site Crédits