Particule chargée dans le champ de deux charges égales

1. Introduction

Ce document traite le mouvement d'une particule chargée dans le champ créé par deux charges ponctuelles identiques et fixées. Le mouvement est considéré dans un plan contenant les deux charges. Il s'agit du problème des trois corps en interaction newtonienne, avec deux corps fixés.

Pour l'intégration numérique des équations différentielles du mouvement, on utilise la fonction scipy.integrate.odeint.

2. Équations du mouvement

On se place dans un plan contenant les deux charges, dans lequel on utilise un repère (0xy). Les deux charges A et B sont placées sur l'axe x respectivement en x=-a et x=+a. Le potentiel électrostatique dans le plan est :

V (x, y) = \frac{q_{0}}{4 π ε_{0}} (\frac{1}{\sqrt{{(x - a)}^{2} + y^{2}}} + \frac{1}{\sqrt{{(x + a)}^{2} + y^{2}}}) (1)

Pour le calcul numérique, on utilise le potentiel réduit suivant, qui suppose q₀>0 :

V (x, y) = \frac{1}{\sqrt{{(x - 1)}^{2} + y^{2}}} + \frac{1}{\sqrt{{(x + 1)}^{2} + y^{2}}} (2)

Les composantes du champ électrique sont :

- \frac{\partial V}{\partial x} = \frac{x - 1}{({(x - 1)}^{2} + y^{2})^{\frac{3}{2}}} + \frac{x + 1}{({(x + 1)}^{2} + y^{2})^{\frac{3}{2}}} (3)

- \frac{\partial V}{\partial y} = \frac{y}{({(x - 1)}^{2} + y^{2})^{\frac{3}{2}}} + \frac{y}{({(x + 1)}^{2} + y^{2})^{\frac{3}{2}}} (4)

Le système d'équations du mouvement d'une particule chargée négative (b>0) est :

\begin{matrix} \frac{d x}{d t} = v_{x} & (5) \\ \frac{d v_{x}}{d t} = b \frac{\partial V}{\partial x} & (6) \\ \frac{d y}{d t} = v_{y} & (7) \\ \frac{d v_{y}}{d t} = b \frac{\partial V}{\partial y} & (8) \end{matrix}

L'intégrale première (énergie mécanique) est :

E = \frac{1}{2} (v_{x}^{2} + v_{y}^{2}) - b V (x, y) (9)

3. Intégration numérique

3.a. Fonction odeint

On utilise la fonction scipy.integrate.odeint. Cette fonction utilise le solveur LSODA de la bibliothèque fortran ODEPACK. La méthode utilisée par ce solveur est une méthode à pas multiple et à ordre variable.

import math
import numpy
from matplotlib.pyplot import *
import scipy.integrate

La fonction suivante définit le système d'équations, avec z[0]=x, z[1]=vx, z[2]=y, z[3]=vz :

def equation(z,t):
    x1 = z[0]-1
    x2 = z[0]+1
    yc = z[2]*z[2]
    d1 = math.pow(x1*x1+yc,1.5)
    d2 = math.pow(x2*x2+yc,1.5)
    return [z[1],-x1/d1-x2/d2,z[3],-z[2]/d1-z[2]/d2]

On choisit le temps maximal du calcul et la période d'échantillonnage puis on crée le tableau des instants :

tmax=50
te=0.01
t = numpy.arange(start=0,stop=tmax,step=te)

La méthode numérique utilisée par odeint effectue une évaluation de l'erreur commise à chaque pas de calcul (passage de l'instant $t$ à $t + Δ t$ ). Le pas $Δ t$ est modifié en permanence de manière à maintenir l'erreur locale inférieure à :

ε_{a} + ε_{r} | z | (10)

où $ε_{a}$ est la tolérance absolue, $ε_{r}$ la tolérance relative et $| z |$ la norme des variables, qui est dans le cas présent :

| y | = \sqrt{x^{2} + y^{2} + v_{x}^{2} + v_{y}^{2}} (11)

Plus les tolérances sont faibles, plus le calcul est précis, mais plus il est long. Les valeurs des tolérances sont déterminées empiriquement. On commence par ces valeurs :

atol=1e-6
rtol=1e-6

La condition initiale choisie est une position proche de l'axe Oy, avec une vitesse nulle :

z0 = [0.2,0,5,0]

L'intégration numérique :

z = scipy.integrate.odeint(equation,z0,t,rtol=rtol,atol=atol)

Le résultat est un tableau numpy dont les colonnes correspondent aux différents instants. Dans le cas présent, ce tableau a 4 lignes. Pour récupérer les listes des coordonnées et des vitesses on extrait les colonnes :

x = z[:,0]
vx = z[:,1]
y = z[:,2]
vy = z[:,3]

Voici le tracé de la trajectoire :

figure(figsize=(6,6))
plot(x,y)
plot([1,-1],[0,0],'r.')
xlabel('x')
ylabel('y')
axis([-6,6,-6,6])

figA.pdf

La particule se rapproche de l'axe sous l'effet des lignes de champ qui convergent vers le centre. Lorsque l'axe Ox contenant les charges est franchi, la particule s'éloigne de l'axe puis est rarentie jusqu'à être attirée vers la charge située B en x=1. Elle subit une diffusion par cette charge puis s'éloigne à nouveau avant de revenir et de subir une deuxième diffusion.

On prolonge la durée du calcul pour voir la suite de la trajectoire :

tmax = 200
t = numpy.arange(start=0,stop=tmax,step=te)
z = scipy.integrate.odeint(equation,z0,t,rtol=rtol,atol=atol)
zt = z.transpose()
x = zt[0]
vx = zt[1]
y = zt[2]
vy = zt[3]
figure(figsize=(6,6))
plot(x,y)
plot([1,-1],[0,0],'r.')
xlabel('x')
ylabel('y')
axis([-6,6,-6,6])

figB.pdf

La particule revient régulièrement vers la charge B au voisinage de laquelle elle subit une diffusion qui dévie fortement sa trajectoire. Ce mouvement est très différent de celui qui serait obtenu avec une seule charge : on aurait alors une diffusion dans un champ coulombien, c'est-à-dire une trajectoire conique.

Pour voir si le mouvement est périodique, on trace la composante sur x sur une plus grande durée :

tmax = 2000
t = numpy.arange(start=0,stop=tmax,step=te)
z = scipy.integrate.odeint(equation,z0,t,rtol=rtol,atol=atol)
zt = z.transpose()
x = zt[0]
vx = zt[1]
y = zt[2]
vy = zt[3]
figure(figsize=(12,4))
plot(t,x)
xlabel('t')
ylabel('x')
axis([0,tmax,-6,6])

figC.pdf

On trace aussi l'énergie mécanique pour vérifier la bonne conduite du calcul numérique :

def energie():
    e = numpy.zeros(t.size)
    for i in range(t.size):
        x1 = x[i]-1
        x2 = x[i]+1
        yc = y[i]*y[i]
        d1 = math.pow(x1*x1+yc,0.5)
        d2 = math.pow(x2*x2+yc,0.5)
        e[i] = 0.5*(vx[i]*vx[i]+vy[i]*vy[i])-1.0/d1-1.0/d2
    return e
figure(figsize=(12,4))
plot(t,energie())
xlabel('t')
ylabel('e')

figD.pdf

On constate une instabilité numérique qui se manifeste vers t=1000 par une brusque variation de l'énergie. À partir de cet instant, le résultat du calcul numérique n'est plus du tout valable.