Discussion :

[Mainot]

Bonjour,

J'essaie de résoudre actuellement l'équation de Laplace avec une configuration particulière, j'ai fait le programme suivant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

Néanmoins, j'ai un problème, j'aimerais maintenant réduire l'écart et le mettre à 10**-3 (while ecart>10**-3 à la place de 5*10**-2), mais le programme n'est sûrement pas optimisé et ne s'arrête jamais.
Je ne vois pas comment faire pour l'améliorer...
Auriez-vous des idées ?
Merci d'avance !

[lg_53]

mais le programme n'est sûrement pas optimisé et ne s'arrête jamais.

Ca à mon avis c'est faux. On voit que l'écart se réduit petit à petit. Certes de plus en plus lentement, mais il se réduit, donc si on le laissait tourné assez longtemps il finirait par s'arêter (assez longtemps ca peut etre des heures, voire des jours ! )

Donc à priori ton schéma numérique qui se résume à cette ligne :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Vnew[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1])

est convergeant. Mais il existe des schémas qui convergent plus rapidemment. On peut faire des moyennes sur plus de points, et on peut aussi faire des moyennes pondérées.

Par exemple si tu essaies avec :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Vnew[i,j] = 0.2*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1]+ V[i,j])

ou

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Vnew[i,j] = 0.125*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1] + 4*V[i,j])

ou bien encore

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
Vnew[i,j] = 1/24*(2*V[i-1,j] + 2*V[i+1,j] + 2*V[i,j-1] + 2*V[i,j+1] + 12*V[i,j]
             + V[i-1,j-1] + V[i-1,j+1] + V[i+1,j-1] + V[i+1, j+1] )

Là tu verras qu'on atteint également difficilement 1e-3, mais que ca descend quand même plus vite.

Ensuite, ce morceau de code m'interpèle :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
  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]

et surtout la ligne avec le if. Un schéma itératif ne peut fonctionner que si à chaque fois on procède sur la nouvelle valeur. Là si la condition n'est plus vérifiée sur une partie de ton domaine de calcul, alors ton schéma stagne à cet endroit là. Il ne peut plus gagner en convergence donc. Et si jamais on n'arrive au point au rien n'ai mis à jour, là le schéma fait du sur place ...

Si j'enlève le "if" et que je mets un schéma plus élaboré, j'obtiens une convergence à 1e-3 au bout d'1 ou 2 minutes d'éxécution.

Après peut-etre que cette condition est importante par rapport au sens qu'elle a. Ca c'est un autre problème. Mais pour le schéma numérique on ne peut pas faire ça.

[Mainot]

Merci pour votre réponse et les alternatives proposées !

Concernant le if, au début je ne l'avais pas mis mais les résultats que j'obtenais étaient faux mais je comprends pourquoi cela peut stagner du coup.

Si j'enlève le "if" et que je mets un schéma plus élaboré, j'obtiens une convergence à 1e-3 au bout d'1 ou 2 minutes d'éxécution

Serait-ce possible savoir ce que vous avez mis à la place du if ?

[lg_53]

J'ai rien mis du tout à la place.

Qu'est ce qui ne va pas dans les résultats si tu ne le mets pas ?

[Mainot]

Au début je n'avais pas tout ce bloc ci :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
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]

j'avais juste mis après le

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

print ecart

:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Vnew = V.copy()

.
Mais les résultats que j'obtenais n'étaient pas bons (l'allure de la courbe et les valeurs du potentiel).

En enlevant toute la condition if, pour 10**-3, c'est toujours aussi long chez moi...
Je ne sais pas si c'est normal, à moins que j'aie mal compris ce que vous vouliez dire.

[lg_53]

Normal que c'était faut au départ !

Il ne faut pas écrire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
print ecart
Vnew = V.copy()

mais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
print ecart
V = Vnew.copy()

On veut que V prenent la nouvelle valeur, pour continuer à travailler avec cette nouvelle valeur.

Et je n'ai pas di que j'avais enlevé toute la boucle, mais juste le if :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
for i in range(N): 
        for j in range(N):
               V[i,j] = Vnew[i,j]

ce qui revient donc à faire

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

V = Vnew.copy()