1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
| from __future__ import division
from pylab import *
from scipy import *
from numpy import *
from matplotlib import *
N1=50 # nombre de point selon x et y
N=2*N1+1 # nombre de points total
# les grandeurs sont données en cm
xc1=4
xc2=9
yc1=0
yc2=0
R1=1.75
R2=9
ecart = 1
a, b = linspace(-1, 19, N), linspace(-10, 10, N)
xa, ya = meshgrid(a, b)
V = zeros_like(xa)
for i in range(N):
for j in range(N):
x, y = xa[i,j], ya[i,j]
if (((x-xc1)**2/(R1**2))+((y-yc1)**2/(R1**2)))<=1 : # definit le potentiel dans le conducteur central
V[i,j] = 30
if (((x-xc2)**2/(R2**2))+((y-yc2)**2/(R2**2)))>=1 : # definit le potentiel dans le conducteur exterieur
V[i,j] =0
# permet de tracer le potentiel selon X le long de l'axe de symétrie.
Vnew = V.copy()
while ecart > 5*10**-2:
for i in range(1,N-1):
for j in range(1,N-1):
Vnew[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1])
# critère de convergence
ecart = np.max(np.abs(V - (Vnew))/np.max(V))
print(ecart)
# sauvegarde dans la grille V de la grille calculée
for i in range(N):
for j in range(N):
x, y = xa[i,j], ya[i,j]
if (((x-xc1)**2/(R1**2))+((y-yc1)**2/(R1**2))) > 1 and (((x-xc2)**2/(R2**2))+((y-yc2)**2/(R2**2))) < 1 :
V[i,j] = Vnew[i,j]
figure |
Partager