CCP Modélisation de systèmes physiques ou chimiques PC 2017

Thème de l'épreuve Autour de l'équation de Poisson
Principaux outils utilisés électrostatique, résolution numérique, programmation, équations aux dérivées partielles, équations différentielles, théorème de Gauss, développement limité, méthode d'Euler
Mots clefs équation de Poisson, méthode de Jacobi, méthode de Gauss-Seidel, numpy, condensateur plan, tube cathodique, conditions de Dirichlet

Corrigé

(c'est payant, sauf le début): - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Extrait gratuit du corrigé

(télécharger le PDF)
           

Énoncé complet

(télécharger le PDF)
                                            

Rapport du jury

(télécharger le PDF)
        

Énoncé obtenu par reconnaissance optique des caractères


SESSION 2017

PCMS006

!

!
!

EPREUVE SPECIFIQUE - FILIERE PC!
!!!!!!!!!!!!!!!!!!!!"
!

MODELISATION DE SYSTEMES PHYSIQUES OU CHIMIQUES
Jeudi 4 mai : 8 h - 12 h!
!!!!!!!!!!!!!!!!!!!!"
N.B. : le candidat attachera la plus grande importance à la clarté, à la 
précision et à la concision de la
!"#$%&'()*+,'+-)+%$)#'#$&+./&+$0.)"+1+!.2"!.!+%.+3-'+2.-&+4-'+/.054.!+6&!.+-).+.!!.-!+#7")()%"8+'4+4.+
signalera sur sa copie et devra poursuivre sa composition en expliquant les 
raisons des in'&'$&'9./+3-7'4+
a été amené à prendre.!

!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
!
!
!
!
!
!
!
Les calculatrices sont autorisées
!
!
!
!
!
!
!
Le sujet est composé de deux parties, largement indépendantes.
!
!
!
!
!
!
!
!
!
!
!
!
!
!

1/15

!

Autour de l'équation de Poisson
Ce problème s'intéresse à la résolution numérique de quelques problèmes 
d'électrostatique. Il se
compose de deux parties.
I. Étude de l'équation de Poisson et de différentes méthodes de résolution 
numérique.
II. Deux études de cas : fil infini chargé et mouvement d'un électron entre les 
plaques d'un condensateur.
Les différentes parties sont largement indépendantes.
Un aide-mémoire numpy/matplotlib/pyplot est présent à la fin du sujet.

Partie I - Équation de Poisson
I.1 - Établissement de l'équation
Q1. Rappeler l'équation de Maxwell-Gauss ainsi que la relation entre le champ E 
et le potentiel
électrostatique V. En déduire l'équation de Poisson :
V +

=0.
0

Préciser les noms et les unités usuelles de  et 0 .
Q2. Citer plusieurs situations physiques en dehors de l'électrostatique pour 
lesquelles il existe une
équation analogue.
I.2 - Équation adimensionnée pour un problème plan
On veut résoudre l'équation de Poisson dans une portion de plan P carrée de 
côté L. On pose :
X = x/L, Y = y/L .
Q3. Montrer qu'on peut écrire l'équation sous la forme suivante :
2 V(X, Y) 2 V(X, Y)
+
+  (X, Y) = 0
X 2
Y 2
où  (X, Y) sera exprimé en fonction de , L et 0 .
I.3 - Discrétisation
Afin de résoudre numériquement l'équation de Poisson, on va utiliser un 
maillage de P, de pas
h = 1/N, et on va transformer les dérivées partielles par des différences entre 
les valeurs de V aux
différents points du maillage (on parle aussi des noeuds du maillage). La 
figure 1 (page suivante)
représente le maillage de P pour N = 5.

2/15

ZE
î

X,- 1 X

Figure 1 -- Maillage de 73 pour N = 5

Q4. En faisant un développement limité à l'ordre 2 autour du point de 
coordonnées (X,, Yj), montrer

' n t rim rl l rde ôZV + ôZV en ce oint so s la forme s i ante'
quo peu exp e avaeu ôX2 ôY2 p u uv .
âZv â2V V(Xi + h, Yj) + V(Xj _ h, Yj) + V(Xï, Yj + h) + V(Xj, Yj _ h) _ 4V(Xj, 
Yj) +O(h)
_ +_ = _ .
ôX2 ôY2 h2

QS. Comme X, : ih et Yj : jh, on note désormais V(i, j) le potentiel V(X,, Y ]) 
en un point (X,, Y ])
du maillage. Montrer alors qu'on peut écrire l'équation de Poisson sous la 
forme suivante :

V(i + l,j) + V(i -- l,j) + V(i,j+ l) + V(i,j-- l) -- 4V(i, j) +p"(i, j) = 0 (l)
p"(i, j) étant une fonction à définir en fonction de p, L, 80 et h.

I.4 - Résolution

La fonction p"(i , j) étant connue, on montre en mathématiques que la solution 
de l'équation de Pois--
son est unique si on fixe les conditions aux limites sur la frontière 77 du 
domaine 50. Ces conditions
sont essentiellement de deux types :
-- on impose le potentiel en tout point de ?" (conditions de Dirichlet),
-- on impose une condition sur les dérivées partielles de V en tout point de 9" 
(conditions de
Neumann).
Dans ce problème, on ne va considérer que des conditions de Dirichlet.

La frontière ? contient naturellement les points du bord de 50 (donc 
appartenant aux quatre côtés du
carré), mais elle peut aussi contenir certains points à l'intérieur de 73 où le 
potentiel est fixé en raison
de la présence d'électrodes.

L'ensemble des points de coordonnées (i, j) est donc composé de deux 
sous-ensembles :
-- ceux dont le potentiel est connu, appartenant à la frontière 7" ,
-- ceux dont le potentiel est inconnu, appartenant à 79 mais pas à 77 (donc 
dans PW").

3/15

Méthode de Jacobi
À partir de l'équation (1), on peut exprimer :
1
V(i, j) = (V(i + 1, j) + V(i - 1, j) + V(i, j + 1) + V(i, j - 1) +  (i, j)) .
4

(2)

La résolution s'effectue alors en deux étapes.
-- Initialisation
a) On fixe le potentiel des points de F à la valeur imposée physiquement (bords 
et électrodes).
b) On donne aux points de potentiel inconnu, donc appartenant à P\F , une 
valeur arbitraire
V0 (i, j), en général nulle.
-- Itérations
On calcule une nouvelle valeur V1 (i, j) des potentiels en appliquant 
l'équation (2) pour tous
les points de P\F , tandis que V1 (i, j) = V0 (i, j) pour les points de F .
Le processus est répété jusqu'à obtenir des valeurs du potentiel quasiment 
stables. En notant
k le nombre d'itérations, on a donc pour le point de coordonnées (i, j) 
n'appartenant pas à la
frontière :
1
Vk+1 (i, j) = (Vk (i + 1, j) + Vk (i - 1, j) + Vk (i, j + 1) + Vk (i, j - 1) +  
(i, j)) .
(3)
4
La convergence de la méthode est vérifiée à l'aide du critère de convergence ek 
, défini par :

1 
ek =
(Vk+1 (i, j) - Vk (i, j))2 .
(4)
N 2 i, j
Le calcul sera stoppé au bout de k itérations, quand ek deviendra inférieur à 
un seuil de convergence
 fixé arbitrairement.
Implémentation informatique
On va utiliser la bibliothèque numpy permettant une utilisation simple des 
tableaux de flottants à deux
dimensions ; un aide-mémoire est disponible en fin de sujet.
Le chargement des bibliothèques classiques est assuré par les lignes suivantes :
# importation des biblioth èques
import numpy as np
import matplotlib . pyplot as plt
import math
On supposera que les tableaux numpy suivants, utilisés comme arguments dans les 
fonctions à définir
dans les questions qui suivent, ont pour signification :
-- V[i,j], (i, j)  0 . . . N2 : tableau courant du potentiel en un point de P,
-- rhos[i,j], (i, j)  0 . . . N2 : tableau contenant la densité de charge  en 
un point de P,
-- frontiere[i,j], (i, j)  0 . . . N2 : tableau de booléens indiquant si le 
point de coordonnées (i, j) appartient ou non à F . En particulier, tous les 
points du bord du domaine seront tels
que frontiere[i,j]==True.

4/15

Q6. Écrire la fonction nouveau_potentiel(V, rhos, frontiere, i, j) retournant 
la nouvelle valeur du potentiel au point (i, j)  0 . . . N2 selon l'équation 
(3).
Q7. Montrer que pour modifier toutes les valeurs contenues dans V[i,j] pendant 
une itération, il
est nécessaire de disposer d'une copie de ce tableau.
On rappelle que l'attribut shape permet de récupérer les dimensions d'un 
tableau numpy.
Q8. Écrire la fonction itere_J(V, rhos, frontiere) modifiant la totalité du 
tableau V[i,j]
lors d'une seule itération et retournant l'erreur calculée conformément à 
l'équation (4).
Q9. Écrire la fonction poisson(f_iter, V, rhos, frontiere, eps) ayant pour 
premier argument une fonction du même type que celle définie à la question 
précédente, pour dernier
argument eps le seuil arbitraire de convergence  et dont le rôle est de 
modifier le tableau des
potentiels V[i,j] jusqu'à convergence.
I.5 - Améliorations
Méthode de Gauss-Seidel
C'est une modification de la méthode de Jacobi, pour laquelle on montre que la 
convergence est
légèrement plus rapide. Supposons que l'on balaye le tableau des potentiels 
selon les indices i et j
croissants : dans ces conditions, les points situés à gauche et en dessous du 
point courant ont déjà
été calculés. On va utiliser ces nouvelles valeurs, probablement plus proches 
de la solution, dans la
formule permettant le calcul de Vk+1 (i, j). Ceci donne l'algorithme de 
Gauss-Seidel :
1
Vk+1 (i, j) = (Vk (i + 1, j) + Vk+1 (i - 1, j) + Vk (i, j + 1) + Vk+1 (i, j - 
1) +  (i, j)) .
4

(5)

Q10. Montrer qu'il n'est plus nécessaire de copier le tableau V[i,j] pour la 
mise à jour lors d'une
itération en utilisant l'équation (5). Faut-il modifier la fonction 
nouveau_potentiel pour
passer de la méthode de Jacobi à celle de Gauss-Seidel ?
Q11. Écrire la fonction itere_GS(V, rhos, frontiere) modifiant la totalité du 
tableau V[i,j]
lors d'une seule itération et retournant l'erreur calculée conformément à 
l'équation (4).
Méthode de Gauss-Seidel adaptative
Les méthodes de Jacobi et de Gauss-Seidel n'utilisent pas la valeur de Vk (i, 
j) pour calculer Vk+1 (i, j).
La méthode de sur-relaxation (Successive Over Relaxation method) consiste à 
calculer la nouvelle
valeur d'un noeud comme une combinaison linéaire de la valeur courante et de 
celle donnée par le
schéma de Gauss-Seidel. En introduisant le paramètre de relaxation , on a alors 
:
Vk+1 (i, j) = (1 - )Vk (i, j) +

(Vk (i + 1, j) + Vk+1 (i - 1, j) + Vk (i, j + 1) + Vk+1 (i, j - 1) +  (i, j)) . 
(6)
4

L'étude mathématique de cette relation permet de montrer les résultats suivants 
:
-- la méthode converge uniquement si 0 <  < 2 et elle converge plus rapidement 
que la méthode
de Gauss-Seidel si 1 <  < 2,
-- il existe une valeur optimale de  qui permet la convergence avec un nombre 
d'itérations en
O(N) pour une valeur de  fixée.

5/15

Pour la résolution de l'équation de Poisson envisagée dans ce problème 
(conditions de Dirichlet sur
un maillage carré), on montre que la valeur optimale opt est :
opt =

2
.
1 + /N

(7)

Q12. Écrire la fonction nouveau_potentiel_SOR(V, rhos, frontiere, i, j, omega) 
retournant la nouvelle valeur du potentiel au point (i, j) selon l'équation (6).
Q13. Écrire la fonction itere_SOR(V, rhos, frontiere) optimale modifiant la 
totalité du tableau V[i,j] lors d'une seule itération et retournant l'erreur 
calculée conformément à l'équation (4).
La résolution du problème peut alors se faire par un appel de la forme
poisson(itere_SOR, V, rhos, frontiere, eps),
les tableaux carrés V, rhos, frontiere étant de dimensions convenables pour 
représenter un
maillage comportant (N + 1)2 noeuds.
Q14. Quelle est la complexité temporelle de l'appel précédent quand  = opt ?
La figure 2 représente, pour  = 10-4 , la durée d'exécution T (en secondes) en 
fonction de N.
Cette courbe est-elle en accord avec la complexité temporelle attendue ?
Quelle serait la durée d'exécution pour N = 1 000 ? Commenter.
!"#$%&''(')*+),(-./0"1"2"#"3"0

Figure 2 ­ Durée d'exécution (en secondes) de poisson(itere_SOR, V, rhos, 
frontiere,
eps) en fonction de N

6/15

I.6 - Détermination du champ électrique
Connaissant le potentiel V[i,j], il est souvent nécessaire de calculer 
numériquement les composantes E x et Ey du champ électrique au niveau des 
noeuds du maillage, qui sont alors conservées dans
deux tableaux Ex et Ey dimensionnés correctement (ce sont donc des tableaux 
carrés de (N + 1)2
éléments).
Q15. Expliquer rapidement comment il serait possible de définir la fonction 
calc_ExEy(Ex, Ey,
V, h), permettant, à partir du tableau V et du pas du maillage h, le 
remplissage des deux
tableaux Ex et Ey passés en arguments.
Remarque : on ne demande pas d'écrire le code de la fonction, juste de décrire 
précisément
les étapes de calcul, ainsi que les différents cas à considérer.

Partie II - Deux études de cas
II.1 - Fil cylindrique chargé uniformément
Étude théorique
On considère dans le vide un fil cylindrique infini d'axe z et de rayon R, 
portant une charge volumique
constante .
Q16. En se plaçant en coordonnées cylindriques d'axe z, montrer par des 
considérations de symétrie
et d'invariance que le champ E = E(r) ur . En déduire la forme des surfaces 
équipotentielles.
Q17. En appliquant le théorème de Gauss, calculer le champ E dans tout 
l'espace. Tracer rapidement
l'allure de E(r) en fonction de r.
Q18. On donne : 0 = 8,85 × 10-12 F.m-1 ,  = 1,00 × 10-5 C.m-3 , R = 5,00 cm. 
Calculer la valeur
maximale de la norme du champ électrique, ainsi que la valeur pour r = 2R.
Étude numérique
Pour pouvoir utiliser la méthode de Gauss-Seidel adaptative, on place le fil 
infini au centre d'une
enceinte de longueur infinie et de section carrée (L × L), portée au potentiel 
nul (figure 3).
Dans la suite, on prendra L = 4R = 20,0 cm.
y
enceinte

L
fil chargé
L
2

2R

V =0
O

L
2

L

x

Figure 3 ­ Fil infini dans une enceinte de section carrée, portée au potentiel 
nul

7/15

Le programme permettant la résolution de ce problème commence ainsi :
# initialisations
eps0 = 8.85e -12
L = 20.0e-2
N = 100 ; h = L/N
rho = 1.00e-5

#
#
#
#

epsilon_0
20 cm
dé finition du maillage
densit é vol. de charge rho

# les tableaux globaux numpy pour le cylindre
rhos_cyl = np.zeros ((N+1,N+1)) # tableau des
V_cyl = np.zeros ((N+1,N+1))
# le potentiel
Ex_cyl = np.zeros ((N+1,N+1)) # la composante
Ey_cyl = np.zeros ((N+1,N+1)) # la composante

charg é
valeurs de rho ''
Ex
Ey

# le tableau dé finissant la fronti ère est initialement
# rempli entiè rement par la valeur False
frontiere_cyl = np.zeros ((N+1,N+1) , bool)
Q19. Écrire la fonction dans_cylindre(x,y,xc,yc,R) retournant un résultat 
booléen indiquant
si le point de coordonnées (x, y) est à l'intérieur ou sur le bord du cercle de 
centre (xc , yc ) et de
rayon R.
Q20. Écrire la fonction initialise_rhos_cylindre(tab_rhos), initialisant le 
tableau
rhos_cyl contenant les valeurs  (i, j) pour les noeuds du maillage.
Q21. Écrire la fonction initialise_frontiere_cylindre(tab_f), mettant à True 
les points
appartenant à la frontière, donc de potentiel fixé.
La résolution numérique avec la méthode de Gauss-Seidel adaptative, utilisant 
les valeurs numériques
précédentes et un seuil de convergence  = 10-5 , mène à la figure 4, où on a 
tracé un réseau de courbes
équipotentielles, le potentiel V et la composante E x du champ le long de l'axe 
de symétrie défini par
y = L/2. En outre, la valeur calculée de Ex_cyl[50, 50] est égale à 
-0.0023321214257521206.
Q22. Commenter le plus complètement possible ces résultats ; on veillera, en 
particulier, à les comparer au modèle théorique (allure des courbes, valeurs 
numériques ...).

 y) · u x en fonction de x pour y = L/2
Figure 4 ­ Équipotentielles, V(x, y) et E(x,
Une autre résolution est effectuée, avec une répartition de charges dans le 
cylindre différente de la
précédente, utilisant la même valeur de la densité volumique  = 10-5 C.m-3 . 
Elle mène aux courbes
de la figure 5 (page suivante).
8/15

 y) · u x en fonction de x pour la nouvelle répartition de
Figure 5 ­ Équipotentielles, V(x, y) et E(x,
charge
Q23. Déduire de ces courbes la répartition de charges dans le cylindre dans 
cette deuxième situation.
Calculer le champ électrique en tout point pour cette répartition de charges 
dans le cas d'un
cylindre infini seul dans l'espace.
On vérifiera que la valeur maximale du champ électrique calculée à l'aide de 
cette modélisation
est compatible avec celle déduite de la figure 5.
II.2 - Mouvement d'un électron dans un tube d'oscilloscope
La figure 6 montre un tube d'oscilloscope de petite dimension, dans lequel des 
électrons émis par la
cathode sont accélérés et déviés vers un écran luminescent. La déviation est 
assurée par le passage des
électrons entre les plaques de deux condensateurs plans : un pour la déviation 
horizontale, l'autre pour
la déviation verticale. L'étude qui suit ne concernera que le condensateur 
responsable de la déviation
verticale.

(!)"*

!"#$%&'

Figure 6 ­ Petit tube d'oscilloscope, de longueur d'environ 20 cm
On modélise la trajectoire d'un électron de la façon suivante (figure 7 page 
suivante, où la zone de
déviation est grisée) :
-- on négligera l'effet de la pesanteur,
-- émis à vitesse nulle par effet thermo-électronique au niveau de la cathode 
portée au potentiel
nul, l'électron est accéléré à l'aide d'une tension V0 > 0 afin d'acquérir à 
l'entrée de la zone
de déviation une vitesse v0 ,
-- pendant son trajet dans la zone de déviation, il est soumis à un champ 
électrique E lié aux
potentiels ±V p des plaques du condensateur, de longueur  et séparées par une 
distance d,
-- poursuivant son mouvement, il arrive sur la surface de l'écran à une 
distance y s de l'axe x,
l'écran étant situé à la distance D du centre du condensateur.
9/15

Figure 7 -- Schéma de la zone de déviation

Valeurs numériques :

masse d'un électron m = 9,11 >< 10_31 kg ; charge élémentaire e = 1,60 X 10_19 C

V0 =950V; Vp =180V; D=7,00cm; d=2,00cm; EUR=4,00cm.
Étude physique

Q24. En appliquant la conservation de l'énergie, calculer la vitesse vo de 
l'électron à l'entrée de la
zone de déviation. Faire l'application numérique. Commenter.

Q25. On modélise les plaques de déviation comme un condensateur sans effets de 
bord : le champ
électrique est donc considéré comme nul si |x| > 6/2 et uniforme si lxl EUR 
5/2, ses lignes de
champ étant parallèles à l'axe Oy. Exprimer le champ Ê entre les plaques en 
fonction de Vp et

d.

Q26. On suppose que la vitesse d'entrée de l'électron dans la zone de déviation 
est 70 : vo üx. En
appliquant les lois de la mécanique, établir l'équation de la trajectoire de 
l'électron entre les
plaques (pour |x| < 6/2).

2EURV EUR
Q27. Montrer que l'équation de la trajectoire pour x > EUR / 2 est donnée par : 
y : dp2 x.
m vo
, . , , ,, Vp EUR D
En déduire que l ordonnee du spot surl ecran est : ys : Î >< Î'
0

Faire l'application numérique pour y,.
Étude numérique

Pour savoir si la modélisation précédente est pertinente, on va envisager une 
détermination numérique
de la trajectoire de l'électron. Pour cela, on place le condensateur de 
déviation dans une enceinte
carrée au potentiel nul de côté L = 10,0 cm, le centre du condensateur étant à 
3,0 cm du bord gauche
de l'enceinte (figure 8 page suivante). Les autres caractéristiques électriques 
et géométriques sont les
mêmes que précédemment.

La résolution numérique se déroule alors en deux étapes :
-- calcul du potentiel et du champ électrique par la méthode de Gauss--Seidel 
adaptative dans
l'enceinte,
-- calcul de la trajectoire de l'électron à l'aide de la méthode d'Euler.

10/15

T \

plaques de déviatidn

Figure 8 -- Modélisation des plaques de déviation dans l'enceinte

Le programme permettant cette résolution numérique commence ainsi :

L = 0.100 ; N = 100 ; h : L/N

V0 = 950 ; Vp : 180

m : 9.11e--31 ; e : 1.60e--19

rhos_osc : np.zeros((N+l,N+l))

V_osc : np.zeros((N+l,N+l))

Ex_osc : np.zeros((N+l,N+l))

Ey_osc : np.zeros((N+l,N+l))

frontiere_osc : np.zeros((N+l,N+l), dtype=bool)

Calcul du potentiel et du champ électrique dans l'enceinte

Q28. Quelles sont les valeurs qui doivent être contenues dans le tableau rhos 
pour le problème

considéré ?

Q29. Écrire la fonction initialise_frontiere_condensateur(tab_V, tab_f), 
permettant
l'initialisation des tableaux V_osc et frontiere_osc à l'aide de la ligne de 
code suivante :

initialise_frontiere_condensateur(V_osc, frontiere_osc)

Pour pouvoir utiliser la méthode d'Euler, il est nécessaire de pouvoir calculer 
les composantes Ex(x, y)
et Ey(x, y) du champ E pour x EUR [O, L[ et y EUR [O, L[. Cependant, la méthode 
de résolution (asso--
ciée à la fonction calc_ExEy définie dans la question 15) ne permet de calculer 
les composantes

Ex_osc [i , j] et Ey_osc [i , j] qu'aux noeuds du maillage.

11/15

Soit un point P de coordonnées (x, y)  [0, L[×[0, L[. Ce point est dans la 
cellule (i, j), où i = x/h et
j = y/h. Posons r x = x - ih et ry = y - jh.
Q30. Montrer alors que :
E x (x, y)  Ex[i,j] +
((Ex[i+1,j] - Ex[i,j]) * rx + (Ex[i,j+1] - Ex[i,j]) * ry) / h .
Écrire de même la formule permettant de calculer Ey (x, y). En déduire qu'il 
est possible de
calculer les composantes du champ en tout point de P.
On supposera dans la suite que les fonctions val_Ex(Ex,Ey,x,y,h) et 
val_Ey(Ex,Ey,x,y,h)
sont définies et retournent les valeurs des composantes du champ électrique 
pour le point M de coordonnées (x, y) calculées à l'aide des formules 
précédentes.
Calcul de la trajectoire par la méthode d'Euler
Q31. Montrer que les équations permettant de décrire le mouvement de l'électron 
par la méthode
d'Euler sont les suivantes, avec t comme petit incrément temporel et x, v x , 
y, vy les variations pendant t des grandeurs x, v x , y, vy :

y = vy t

 x = v x te
e
.

vy = - Ey (x, y)t
 v x = - E x (x, y)t
m
m
Q32. Compte tenu du changement de la position d'origine du repère (imposée par 
la résolution numérique de l'équation de Poisson, figure 8), quelles sont les 
conditions initiales du mouvement
de l'électron ?
Déterminer t pour calculer environ 200 points successifs le long de la 
trajectoire. Faire l'application numérique.
Q33. En tenant compte des réponses aux questions précédentes, compléter le code 
d'initialisation
des variables de la simulation (****** dans le code suivant) :
Npts = 200 # nombre de points pour le tracé de la trajectoire
v0 = ****** # vitesse initiale de l'é lectron
dt = ****** # incrément temporel
# tableaux des coordonn ées x et y de l'é lectron
lx = np.zeros(Npts) ; ly = np.zeros(Npts)
# tableaux des vitesses en x et en y
lvx = np.zeros(Npts) ; lvy = np.zeros(Npts)
# conditions initiales
lx [0] = ****** ; ly [0] = ******
lvx [0] = ******* ; lvy [0] = *******
Q34. Écrire les lignes de code implémentant la boucle de remplissage des 
tableaux lx, ly, lvx,
lvy selon la méthode d'Euler.

12/15

Comparaison théorie/simulation
La figure 9 montre le résultat de la simulation précédente. On y voit le réseau 
de courbes équipotentielles, ainsi que deux trajectoires 1 et 2 , l'une étant 
associée au calcul théorique, l'autre à la
simulation numérique. Chaque trajectoire est constituée de 200 points de calcul 
séparés d'une durée
t.
Q35. Reproduire sommairement sur la copie la figure 9, y ajouter le tracé de 
quelques lignes de
champ orientées dans les différentes parties de la zone de déviation.
Identifier, en le justifiant, chaque trajectoire. Expliquer pourquoi la 
trajectoire 1 est plus
courte que la trajectoire 2 .
À votre avis, peut-on se contenter de l'étude théorique pour prévoir le point 
d'impact de l'électron sur l'écran ? (On attend une réponse chiffrée.)

!

"

Figure 9 ­ Tracé dans la zone de déviation de quelques courbes 
équipotentielles, des trajectoires de
l'électron (théorique et simulée numériquement)

13/15

Aide-mémoire numpy/matplotlib/pyplot

Importation des bibliothèques
Les bibliothèques sont importées de la façon suivante :
import math
import numpy as np
import matplotlib . pyplot as plt
Manipulation des tableaux numpy
La création d'un tableau numpy à deux dimensions dont toutes les valeurs sont 
initialisées à 0 est
faite par l'instruction np.zeros(format), format étant un doublet de la forme 
(n_lignes ,
n_colonnes) :
>> t0=np.zeros ((2 ,3)); print (t0)
[[ 0. 0. 0.]
[ 0. 0. 0.]]
Pour avoir un tableau rempli de 1, on utilise np.ones(format) :
>> t1=np.ones ((2 ,2)); print (t1)
[[ 1. 1.]
[ 1. 1.]]
On peut récupérer le format d'un tableau en demandant son attribut shape, ce 
qui retourne un doublet :
>> print (t0.shape); print (t1.shape)
(2, 3)
(2, 2)
Dans le cas d'un tableau carré, on peut donc récupérer le nombre de lignes, 
égal au nombre de colonnes, en accédant au premier élément du doublet :
>> t1.shape [0]
2
On peut créer un tableau numpy de booléens en ajoutant le type bool. La valeur 
0 est associée à
False, la valeur 1 à True :
>> np.zeros ((2 ,3) , bool)
[[ False False False]
[False False False ]]}
>> np.ones ((2 ,3) ,bool)
[[ True True True]
[ True True True ]]

14/15

L'accès à un élément du tableau a (en lecture ou en modification) se fait par 
a[i,j], les lignes et les
colonnes étant numérotées à partir de 0 :
>>
>>
[[
[

a=np.zeros ((2 ,3)) ; a[0 ,0]=1 ; a[1 ,2]=2
a
1. 0. 0.]
0. 0. 2.]]

Une copie indépendante d'un tableau a se fait à l'aide de np.copy(a)
>>
>>
[[
[
>>
[[
[

b = np.copy(a) ; b[0 ,0]=3 ; b[1 ,1]=5
a
1. 0. 0.]
0. 0. 2.]]
b
3. 0. 0.]
0. 5. 2.]]

FIN

15/15

Extrait du corrigé obtenu par reconnaissance optique des caractères



CCP Modélisation de systèmes physiques
ou chimiques PC 2017 -- Corrigé
Ce corrigé est proposé par Virgile Andreani (ENS Ulm) ; il a été relu par 
JeanJulien Fleck (professeur en CPGE) et Cyril Ravat (professeur en CPGE).

Ce sujet s'intéresse à la résolution numérique de quelques phénomènes 
électrostatiques. Il s'organise en deux parties indépendantes.
· Dans la première, on résout l'équation de Poisson avec des conditions aux 
limites de Dirichlet sur un maillage carré, au moyen de trois variantes de la
méthode de Jacobi. Cette partie se divise à parts égales entre l'établissement
du schéma numérique et son implémentation.
· La deuxième partie propose deux applications de l'algorithme de la première
partie. On étudie d'abord le champ électrostatique généré par un fil chargé. 
Puis
on simule le fonctionnement d'un tube cathodique d'oscilloscope en modélisant
la trajectoire d'un électron qui se déplace dans le champ créé par deux plaques
parallèles.
L'épreuve étant assez longue, il faut veiller à traiter efficacement certaines 
questions calculatoires comme la démonstration des schémas numériques. Elle 
aborde la
résolution des équations différentielles partielles et ordinaires, et utilise 
le résultat des
premières pour simuler les secondes. À part une petite imprécision sur la 
variable h
(qui change de dimension au cours du sujet), l'épreuve est en général claire et 
bien posée, et questionne à plusieurs reprises la validité des résultats 
obtenus. Ces questions
portant sur l'interprétation des résultats sont les plus difficiles car elles 
nécessitent de
l'observation ainsi qu'un bon esprit analytique. Le reste de l'épreuve est de 
difficulté
moyenne, et contient notamment deux études très classiques d'électrostatique et 
de
mécanique.

Indications
Partie I
2 L'exemple classique concerne la gravitation. Pour les autres, il faut 
chercher les
formules faisant intervenir le gradient et la divergence.
4 Commencer par développer V(Xi +
- h, Yj ) et V(Xi , Yj +
- h) à l'ordre 2 en h, puis
les combiner pour atteindre l'expression demandée.
7 Que se passerait-il si l'on modifiait les données du tableau au fur et à 
mesure ?
11 Prendre soin de calculer l'erreur au bon moment : entre le calcul du nouveau
potentiel et l'écrasement de l'ancien.
14 Pour estimer graphiquement la complexité, on peut comparer les temps 
d'exécution pour deux valeurs de N espacées d'un facteur 2.
15 Faire attention à ce qui se passe aux frontières du domaine.
Partie II
20 Attention, la définition de h a changé ici par rapport au début de l'énoncé 
: c'est
maintenant L/N au lieu de 1/N.
22 Chercher les résultats obtenus précédemment dans les graphes. Ne pas oublier 
de
commenter la valeur numérique donnée.
23 Le champ électrostatique est nul au centre, linéaire sur une tranche de 
largeur R/2,
puis hyperbolique.
26 Attention à ne pas se tromper de conditions initiales.
28 Les potentiels des armatures sont fixés.
dEy (ih, jh)
dEx (ih, jh)
et
en fonction
dx
dy
de Ex[i, j], Ex[i+1, j] et Ex[i, j+1].

30 Commencer par établir les expressions de

31 Décrire brièvement ce que réalise la méthode d'Euler.

I. Équation de Poisson

-
1 L'équation de Maxwell-Gauss relie la divergence de E à la densité volumique de
charges électriques  :

-

div E =
0

-
Par ailleurs E est l'opposé du gradient du potentiel électrostatique V :
--
-

E = - grad V
--
De plus, V = div(grad V). En combinant ces trois équations, on obtient bien
V +

=0
0

La densité volumique de charge électrique  s'exprime en C · m-3 . Et l'unité
usuelle de 0 , la permittivité diélectrique du vide, est F · m-1 .
2 On retrouve l'équation de Poisson dans plusieurs domaines de la physique, les

-
grandeurs E , V,  et 0 pouvant respectivement s'identifier à
 -
-

· G = F/m,  = Ep /m,  et -1/(4  G) en gravitation newtonienne. On retrouve

-
4  G + div G = 0,

--
-

G = - grad 

et  = 4  G

-
· j , T, P (puissance de production thermique) et  dans la conduction thermique
en régime stationnaire. On retrouve

-
div j = P,

--
-

j = - grad T

P
=0

et T +

-
· j , n, P (terme de production) et D dans la diffusion de particules, également
en régime stationnaire. On retrouve aussi

-
div j = P,

--
-

j = -D grad n

et n +

P
=0
D

· En prenant la divergence de l'équation d'Euler incompressible en régime per
1 --

manent grad p + div -
v -
v = 0, on obtient aussi une équation de Poisson :

-
-

1
-

p + div U = 0 avec U = div 
v -
v

L'intégralité de cette réponse n'était sûrement pas attendue ! On aurait pu
se contenter des deux ou trois premiers exemples.
3 En développant le laplacien, on obtient
 2 V(x, y)  2 V(x, y) (x, y)
+
+
=0
x2
y 2
0
or, par changement de variable,
 2 V(x, y)
 2 V(X, Y)
=
x2
X2

dX
dx

2

=

 2 V(X, Y) 1
X2
L2

et de façon similaire pour

 2 V(x, y)
, d'où
y 2

 2 V(X, Y)  2 V(X, Y)
+
+  (X, Y) = 0
X2
Y2

avec

 (X, Y) = (x, y)

L2
0

4 On peut exprimer V(Xi +
- h, Yj ) et V(Xi , Yj +
- h) à l'aide d'un développement
limité d'ordre 2 autour de V(Xi , Yj ) :
V(Xi , Yj )
 2 V(Xi , Yj )
h+
X
X2
2
V(X
,
Y
)

V(X
i
j
i , Yj )
V(Xi , Yj +
h+
- h) = V(Xi , Yj ) +
-
Y
Y2
On en déduit que
V(Xi +
- h, Yj ) = V(Xi , Yj ) +
-

h2
+ O(h3 )
2
h2
+ O(h3 )
2

 2 V(Xi , Yj ) 2
h + O(h3 )
X2
 2 V(Xi , Yj ) 2
h + O(h3 )
V(Xi , Yj + h) + V(Xi , Yj - h) = 2 V(Xi , Yj ) +
Y2
En combinant ces équations, on obtient bien
V(Xi + h, Yj ) + V(Xi - h, Yj ) = 2 V(Xi , Yj ) +

V(Xi + h, Yj ) + V(Xi - h, Yj ) + V(Xi , Yj + h) + V(Xi , Yj - h)
2V 2V
+
=
2
2
X
Y
h2
V(Xi , Yj )
-4
+ O(h)
h2
En faisant les calculs à l'ordre supérieur, on se rend compte que les termes
en h3 s'annulent également : on pourrait donc remplacer le O(h) du résultat
final par un O(h2 ).
Pour rappel, f (h) = O(hn ) exprime que f (h) est bornée asymptotiquement par 
hn : il existe un k > 0 tel qu'à partir d'un certain h,
|f (h)| 6 k · hn
5 En insérant le résultat de la question précédente dans celui de la question 
3, on
obtient, avec les notations discrètes de l'énoncé,
V(i + 1, j) + V(i - 1, j) + V(i, j + 1) + V(i, j - 1) - 4 V(i, j)
L2
+
(X
L,
Y
L)
=0
i
j
h2
0
d'où l'on déduit que

 (i, j) = (i h L, j h L)

h 2 L2
0

Profitons-en pour vérifier les unités : h est sans dimension, l'unité de  est
C m2
donc 3
, ce qui se réduit bien à des volts, l'unité de V. Attention, la
m F/m
définition de h, ainsi que sa dimension, change en cours de sujet. Dans la
dernière partie, h est en effet défini par L/N et est donc bien une longueur.
6 Si le point de coordonnées (i, j) appartient à F , ses valeurs à l'itération 
k et k + 1
sont identiques. Dans ce cas, on retourne donc sa valeur précédente. Sinon, on 
calcule
sa nouvelle valeur à l'aide de l'équation (3).