import copy

import numpy as np
import [Link] as plt

# Échelle de discrétisation
N
= 10
format = (N+1, N+1)

"""
Choix d'unité :
* Dans la partie
théorique, le sujet considère que les indices 0 à N (inclus) correspondent à
des longueurs
comprises entre 0 et 1 (on travaille avec les variables d'espace X et Y qui
sont sans
dimension).
* Dans les applications, changement de décor, on travaille avec les variables
d'espaces x et y
(des longueurs) qui varient entre 0 et L.

* J'adopte dans tout le

code la convention des applications. En particulier, le pas de la
discrétisation est défini
par h = L/N.

Incidences :
* tableau 'rhos' (charge volumique) : rhos[i,j] =
h**2*rho(i*h, j*h)/epsilon_0
* tableau 'Ex' et 'Ey' (coordonnées du champ
électrostatique) :
Ex[i,j] = (V[i,j]-V[i+1,j])/h
Ey[i,j] =
(V[i,j]-V[i,j+1])/h
"""

## Estimation de la durée d'exécution

# from
[Link] import linregress as linregress
# lx = [Link]([Link](50, 150, 10))
# y = [0.8,
1.3, 2, 3, 4.5, 5.8, 8, 10.2, 13, 16]
# [Link](lx, [Link](y), 'rv')
# [Link]()
#

# a, b, r, p_v, sigma = linregress(lx, [Link](y))

# [Link](lx, a*lx+b, 'b')

##
Comparaison des trois méthodes sur un exemple arbitraire

V = [Link](format)
rhos =
[Link](format) # cas d'une absence de charge volumique

# Calcul du potentiel par

itérations successives
def nouveau_potentiel(V, rhos, frontiere, i, j):
if
frontiere[i,j]:
valeur = V[i,j] # potentiel fixé
else:
valeur =
(V[i+1,j]+V[i-1,j]+V[i,j-1]+V[i,j+1]+rhos[i,j])/4
return valeur

# Variante avec Cte de

relation (omega)
def nouveau_potentiel_SOR(V, rhos, frontiere, i, j, omega):
return
(1-omega)*V[i,j]+omega*nouveau_potentiel(V, rhos, frontiere, i, j)

# Mesure d'écart
entre deux approximations successives
def erreur(V, V_0):
nb_lig, nb_col = [Link]
N =
nb_lig - 1
ecart_quad = [Link]([Link]((V - V_0)**2))
ecart_quad_moy = ecart_quad/N

return ecart_quad_moy

# Première méthode : schéma de Jacobi

def itere_J(V, rhos,
frontiere):
format = [Link]
V_0 = [Link]()
nb_lig, nb_col = format
for i in
range(nb_lig):
for j in range(nb_col):
V[i,j] = nouveau_potentiel(V_0,
rhos, frontiere, i, j)
err = erreur(V, V_0)
return err

# Deuxième méthode : schéma

de Gauss-Seidel
def itere_GS(V, rhos, frontiere):
nb_lig, nb_col = [Link]
V_0 =
[Link]()
for i in range(nb_lig):
for j in range(nb_col):
V[i,j] =
nouveau_potentiel(V, rhos, frontiere, i, j)
err = erreur(V, V_0)
return err

#
Troisième méthode : schéma de Gauss-Seidel adaptatif
def itere_SOR(V, rhos, frontiere):

nb_lig, nb_col = [Link]

N = nb_lig - 1
omega = 2/(1+[Link]/N)
V_0 = [Link]()

for i in range(nb_lig):
for j in range(nb_col):
V[i,j] =
nouveau_potentiel_SOR(V, rhos, frontiere, i, j, omega)
err = erreur(V, V_0)
return
err

# Résolution approchée de l'équation de Poisson

def poisson(f_iter, V, rhos,
frontiere, eps):
i = 0 # compteur du nb d'itérations
err = 2*eps
while
(err>eps):
i += 1
err = f_iter(V, rhos, frontiere)
print(i)

def
init_potentiel(N):
V = [Link]((N+1,N+1))
# Valeurs imposées du potentiel
for i
in range(N+1):
V[i,0] = i/N
V[0,i] = i/N
V[i,N] = 1 # 1-i/N
V[N,i] = 1 # 1-i/N
return V

def comparatif(N, eps):

format = (N+1, N+1)
# Pas de
charge volumique !
rhos = [Link](format)
# Définition de la frontière (= là où
le potentiel est imposé)
frontiere = [Link](format, dtype=[Link])
for i in
range(N+1):
frontiere[i,0]=True
frontiere[i,N]=True

frontiere[0,i]=True
frontiere[N,i]=True

V = init_potentiel(N)
poisson(itere_J,
V, rhos, frontiere, eps)
[Link](V, origin="lower")

plt.tick_params(labeltop=False, labelbottom=True)

V = init_potentiel(N)

poisson(itere_GS, V, rhos, frontiere, eps)

[Link](V, origin="lower")

plt.tick_params(labeltop=False, labelbottom=True)

V = init_potentiel(N)

poisson(itere_SOR, V, rhos, frontiere, eps)

[Link](V, origin="lower")

plt.tick_params(labeltop=False, labelbottom=True)

# N = 50
# eps = 0.001 - Nb
d'itérations : 187 / 230 / 93
# eps = 0.0005 - Nb d'itérations : 452 / 403 / 104
#
eps = 0.0001 - Nb d'itérations : 1260 / 811 / 117
# eps = 0.00005 - Nb
d'itérations : 1611 / 987 / 122

## Calcul du champ en fonction du potentiel

def
calcul_ExEy(Ex, Ey, V, h):
# cas général
Ex[:-1,:] = (V[:-1]-V[1:])/h
Ey[:,:-1] =
(V[:,:-1]-V[:,1:])/h
# cas particuliers
Ex[-1,:] = Ex[-2,:] # dernière ligne =
avant-dernière ligne
Ey[:,-1] = Ey[:,-2] # dernière colonne = avant-dernière colonne

#
# essai
# N = 5
# format = (N+1,N+1)
# V = init_potentiel(N)
# frontiere = [Link](format,
dtype=[Link])
# for i in range(N+1):
# frontiere[i,0]=True
# frontiere[i,N]=True
#
frontiere[0,i]=True
# frontiere[N,i]=True
# poisson(itere_SOR, V, rhos, frontiere, 0.001)
#
Ex, Ey = [Link](format), [Link](format)
# calcul_ExEy(Ex, Ey, V, h)

## Application 1 :
fil cylindrique dans une enceinte carrée

## Q19
def dans_cylindre(x, y, xc, yc, R):
#
cylindre plein uniformément chargé
return (x-xc)**2+(y-yc)**2<=R**2

def
dans_cylindre_creux(x, y, xc, yc, R):
# variante avec un cylindre creux
d2 =
(x-xc)**2+(y-yc)**2
return (d2<R**2) and (d2>0.25*R**2)

## Q20
def
initialise_rhos_cylindre(tab_rhos):
valeur = rho*h**2/eps0
for i in range(N+1):

for j in range(N+1):
if dans_cylindre(i*h, j*h, L/2, L/2, L/4):

tab_rhos[i,j] = valeur

def initialise_rhos_cylindre(tab_rhos, condition):

# variante perso
valeur = rho*h**2/eps0
for i in range(N+1):
for j in
range(N+1):
if condition(i*h, j*h, L/2, L/2, L/4):
tab_rhos[i,j] =
valeur

## Q21
def initialise_frontiere_cylindre(tab_f):
for i in range(N+1):

tab_f[i,0]=True
tab_f[i,N]=True
tab_f[0,i]=True
tab_f[N,i]=True

## Mise en œuvre
eps0 = 8.85e-12
L = 20.0e-2
N = 100
h = L/N
rho = 1.00e-5

def
mise_en_oeuvre(condition):
format = (N+1,N+1)
rhos_cyl = [Link](format)
V_cyl =
[Link](format)
Ex_cyl, Ey_cyl = [Link](format), [Link](format)
frontiere_cyl =
[Link](format, [Link])
# Initialisation charges et potentiel

initialise_rhos_cylindre(rhos_cyl, condition)
initialise_frontiere_cylindre(frontiere_cyl)
 # (potentiel nul sur les points de la frontière, donc
# inutile d'initialiser le
tableau V_cyl !)
# Calcul approché du potentiel
poisson(itere_SOR, V_cyl, rhos_cyl,
frontiere_cyl, 0.00001)
# Calcul approché du champ électrostatique

calcul_ExEy(Ex_cyl, Ey_cyl, V_cyl, h)

# Trois graphes sur une même figure
fig, (ax1,
ax2, ax3) = [Link](1, 3, figsize=(9, 3))
# Lignes de niveau du potentiel
x, y =
[Link]([Link](0, L, N+1), [Link](0, L, N+1))
[Link](x, y, V_cyl, 25,
colors='grey')
ax1.set_title("Équipotentielles", y=1.05)

ax1.set_aspect('equal')
ax1.set_xlim(0, L)
ax1.set_ylim(0, L)

positions_etiquettes = [Link](0, L+0.01, 0.05)

ax1.set_xticks(positions_etiquettes)

ax1.set_yticks(positions_etiquettes)
etiquettes = ["{:d}".format(i) for i in
range(0, 21, 5)]
ax1.set_xticklabels(etiquettes)
ax1.set_yticklabels(etiquettes)

ax1.set_xlabel("x (cm)")
ax1.set_ylabel("y (cm)")
# Potentiel à
mi-hauteur
[Link]([Link](0, L, N+1), V_cyl[:,N//2])
[Link]()

ax2.set_title("V(x, y=L/2) en kV", y=1.05)

ax2.set_xticklabels(etiquettes)

ax2.set_xlabel("x (cm)")
ax2.set_yticklabels(["{:g}".format(0.1*i) for
i in range(0, 19, 2)])
ax2.set_xlim(0, L)
# Composante en x du champ à mi-hauteur

[Link]([Link](0, L, N+1), Ex_cyl[:,N//2])

[Link]()
ax3.set_title("Ex(x,
y=L/2) en kV/m", y=1.05)
ax3.set_xticklabels(etiquettes)
ax3.set_xlabel("x
(cm)")
ax3.set_yticklabels(["{:d}".format(10*i) for i in range(-3, 4)])

ax3.set_xlim(0, L)
plt.tight_layout()

# Patience, l'exécution est assez longue !

mise_en_oeuvre(dans_cylindre)
mise_en_oeuvre(dans_cylindre_creux)

## Application 2 : tube
d'oscilloscope

L = 0.1
N = 100
h = L/N
m = 9.11e-31
e = 1.6e-19
V0, Vp = 950, 180
D =
7e-2
d = 2e-2
ell = 4e-2

format = (N+1,N+1)

V_osc = [Link](format)
Ex_osc, Ey_osc =
[Link](format), [Link](format)

## Q28
rhos_osc = [Link](format)
# et on ne changera rien
: il n'y a pas de charge volumique ici

## Q29
frontiere_osc = [Link](format,
[Link])
def initialise_frontiere_condensateur(tab_V, tab_f):
# la frontière est
constituée du bord du tableau
# où le potentiel est déjà fixé à 0
for i in
range(N+1):
tab_f[i,0] = True
tab_f[i,N] = True
tab_f[0,i] = True

tab_f[N,i] = True
# et aussi des plaques du condensateur
# où le potentiel est fixé
à ±Vp
for i in range(10, 50):
tab_f[i,40] = True
tab_V[i,40] = -Vp

tab_f[i,60] = True
tab_V[i,60] = Vp

initialise_frontiere_condensateur(V_osc,
frontiere_osc)
poisson(itere_SOR, V_osc, rhos_osc, frontiere_osc, 0.00001)
calcul_ExEy(Ex_osc,
Ey_osc, V_osc, h)

# lignes de niveau du potentiel

"""
Problème qui se posait
déjà mais qui passait inaperçu à cause de la symétrie
dans l'application précédente
: il faut remettre les coordonnées dans l'ordre.

On a définit le tableau V et tous les

tableaux analogues en prenant le
premier indice pour représenter les abscisses x et le second
indice pour
représenter les ordonnées y.

Pbm ! Pour une matrice, le premier indice indique

le numéro de ligne et
le second indice le numéro de colonne. Autrement dit, les abscisses
sont
portées sur l'axe vertical et les ordonnées sur l'axe horizontal...

Conséquence : nous allons transposer le tableau V.

"""

def
equipotentielles():
x, y = [Link]([Link](0, L, N+1), [Link](0, L, N+1))

import matplotlib as mpl

[Link]['contour.negative_linestyle']='solid'
 fig, ax =
[Link]()
ax.set_title("Équipotentielles et lignes de champ")
#
potentiels positifs
[Link](x, y, V_osc.transpose(), [Link](Vp/15, Vp, 15),
colors='blue')
# potentiels négatifs
[Link](x, y, V_osc.transpose(),
[Link](-Vp, -Vp/15, 15), colors='red')
# potentiel nul

[Link]([0,0.10],[.05,.05],'k')
xy_ticks = [0, 0.02, 0.04, 0.06, 0.08, 0.1]

xy_etiquettes = ["{:d}".format(i) for i in range(0, 11, 2)]

ax.set_xlabel("x
(cm)")
ax.set_ylabel("y (xm)")
ax.set_xticks(xy_ticks)

ax.set_yticks(xy_ticks)
ax.set_xticklabels(xy_etiquettes)

ax.set_yticklabels(xy_etiquettes)
ax.set_aspect('equal')
ax.set_xlim(0,
0.10)
ax.set_ylim(0, 0.10)
[Link]()
return ax

ax =
equipotentielles()
[Link](x, y, Ex_osc.transpose(), Ey_osc.transpose(),
color='green', density=0.75)

fig, ax = [Link]()
ax.set_title("Intensité
du champ électrostatique (kV/m)")
norme_E_osc = [Link](Ex_osc**2+Ey_osc**2)
xy_ticks =
[0, 0.02, 0.04, 0.06, 0.08, 0.1]
xy_etiquettes = ["{:d}".format(i) for i in range(0,
11, 2)]
ax.set_xlabel("x (cm)")
ax.set_ylabel("y
(xm)")
ax.set_xticks(xy_ticks)
ax.set_yticks(xy_ticks)
ax.set_xticklabels(xy_etiquettes)
a
x.set_yticklabels(xy_etiquettes)
cax = [Link](x, y, norme_E_osc.transpose(),
[Link](0, 2.5e4, 30))
cbar = [Link](cax, ticks=[i*5000 for i in
range(6)])
[Link].set_yticklabels(["0", "5", "10",
"15", "20", "25"])
[Link]('scaled')
[Link](0,
0.10)
[Link](0, 0.10)

fig, ax = [Link]()
ax.set_title("Un champ uniforme ?
(kV/m)")
norme_E_osc = [Link](Ex_osc**2+Ey_osc**2)
xy_ticks = [0, 0.02, 0.04, 0.06, 0.08,
0.1]
xy_etiquettes = ["{:d}".format(i) for i in range(0, 11,
2)]
ax.set_xlabel("x (cm)")
ax.set_ylabel("y
(xm)")
ax.set_xticks(xy_ticks)
ax.set_yticks(xy_ticks)
ax.set_xticklabels(xy_etiquettes)
a
x.set_yticklabels(xy_etiquettes)
cax = [Link](x, y, norme_E_osc.transpose(),
[Link](1.6e4, 1.905e4, 30))
cbar = [Link](cax, ticks=[16000, 16500, 17000, 17500,
18000, 18500, 19000])
[Link].set_yticklabels(["16", "16,5",
"17", "17,5", "18", "18,5",
"19"])
[Link]('scaled')
[Link](0, 0.06)
[Link](0.035, 0.065)

##
Q30
def valeur_champ(E, x, y, h):
# Calcul du champ en tout point de l'espace par
interpolation affine
i, j = x//h, y//h
rx, ry = x-i*h, y-j*h
dEx, dEy =
E[i+1,j]-E[i,j], E[i,j+1]-E[i,j]
return E[i,j]+(dEx*rx + dEy*ry)/h

def val_Ex(Ex, Ey, x,

y, h):
return valeur_champ(Ex, x, y, h)

def val_Ey(Ex, Ey, x, y, h):

return
valeur_champ(Ey, x, y, h)

## Q33
Npts = 200 # nb de points calculés pour le
tracé de la trajectoire
v0 = [Link](2*e*V0/m) # vitesse initiale de l'électron
dt =
L/(Npts*v0) # pas temporel pour le schéma d'Euler

## Q34
def
trajectoire_theorique(x):
K = e*Vp/(m*d*v0**2)
# Entre les armatures, la charge est
soumise à un champ uniforme
avec_force = K*(x+ell/2)**2*(-ell/2<x)*(x<ell/2)
#
En dehors du condensateur, la charge suit un mouvement rectiligne uniforme
sans_force =
2*K*ell*x*(x>ell/2)
return 0.05+avec_force+sans_force

def schema_Euler(Ex, Ey, h):

# Coordonnées de l'électron
x, y = [Link](Npts), [Link](Npts)
# Projetés de
la vitesse de l'électron
vx, vy = [Link](Npts), [Link](Npts)
# Condition
initiales
y[0] = 5e-2
vx[0] = v0
# x[0] et vy[0] sont déjà égaux à 0 !
for
k in range(1, Npts):
x[k] = x[k-1] + vx[k-1]*dt
y[k] = y[k-1] + vy[k-1]*dt

vx[k] = vx[k-1] - e/m*val_Ex(Ex, Ey, x[k-1], y[k-1], h)*dt

 vy[k] = vy[k-1] -
e/m*val_Ey(Ex, Ey, x[k-1], y[k-1], h)*dt
return x, y, vx, vy

ax =
equipotentielles()
x_num, y_num, vx_num, vy_num = schema_Euler(Ex_osc, Ey_osc,
h)
[Link](x_num, y_num, "blue")
x_th = x[0]
y_th =
trajectoire_theorique(x_th+D-L)
[Link](x_th, y_th, "green")

Powered by TCPDF ([Link])