Programme trop long

**Mainot** · 09/12/2017, 11h57

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** · 09/12/2017, 19h22

Envoyé par Mainot

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** · 09/12/2017, 20h01

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** · 09/12/2017, 22h37

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** · 09/12/2017, 23h06

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** · 10/12/2017, 11h01

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()

**Mainot** · 10/12/2017, 14h26

Ahh je vois !

J'ai essayé du coup mais, les résultats que j'obtiens (d'un point de vu physique) ne sont toujours pas bons...
Je pense donc laisser la condition if et mettre 5*10**-2.

Merci encore pour l'aide, si vous avez une autre idée, n'hésitez pas !

**lg_53** · 10/12/2017, 15h40

J'aurais peut etre du commencé par là : tu cherches à représenter quoi physiquement, et quel est le résultat que t'attends à obtenir ?

**Mainot** · 10/12/2017, 15h54

J'essaye de représenter le potentiel (en volts) le long de l'axe x dans cette configuration (les zones grisées sont des conducteurs, le conducteur central est à +15V et le conducteur externe est à -15V).
Nom : Capture d’e?cran 2017-12-10 a? 15.47.37.png
Affichages : 496
Taille : 111,0 Ko

Nom : Capture d’e?cran 2017-12-10 a? 15.47.37.png
Affichages : 496
Taille : 111,0 Ko

Puis de comparer mes résultats expérimentaux avec les résultats théoriques que me donnent le programme :
Nom : TPbonus.png
Affichages : 502
Taille : 338,7 Ko

Nom : TPbonus.png
Affichages : 502
Taille : 338,7 Ko

(précision à 5*10**-2)

Avec une précision plus élevée, normalement, les résultats théoriques et expérimentaux "collent" plus (notamment sur la partie droite de la courbe).

**lg_53** · 10/12/2017, 18h47

Plusieurs remarques qui m'ont conduit au code en bas de mo message que tu essaieras :

1) Evite d'importer des packages avec "import *" Si tu as 40 packages, et que dans le codes il y a une fonction que tu connais pas, et bien on doit chercher dans les 40 packages pour savoir d'ou elle vient. Sans compter que ça peut engendrer des conflits si deux packages définissent chacun une fonction qui porterait le même nom. Bref, à proscrire.

2) Il faut factoriser ton code lorsque tu le peux. Là tu as plusieurs fois une même ligne qui apparait, donc tu peux la factoriser en la mettant dans une fonction. Ca rend le code plus lisible, et plus maintenable (sans compter qu'en plus la fonction peut porter un nom claire qui va ajouter à la sémantique du problème que tu traites)

3) Pour en revenir plus au problème, je comprends maintenant pourquoi tu mets une condition if comme ça dans ta boucle. Sauf que du coup c'est pas très malin de calculer Vnew sur tous le domaine, et dire V=Vnew seulement là où je n'impose pas la valeur (les valeurs imposés étant 0 et 30 dans chacun des composants). C'est plus malin de ne calculer Vnew que si tu es dans la partie où tu as qqch à calculer !

Avec ces 3 remarques j'abouti au code ci dessous, qui s'approche bien du résultat que tu veux. Précision 1e-4, atteinte en moins de 1600 itérations (compter 2-3 minutes d'éxécution). Si je prend un schéma plus élaboré, le résultat voulu est atteint au bout de moins d'itérations, mais pas sûre que ça gagne beaucoup de temps.

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
import numpy as np
from matplotlib import pyplot as plt
 
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
 
a, b = np.linspace(-1, 19, N), np.linspace(-10, 10, N)
xa, ya = np.meshgrid(a, b) 
V = np.zeros_like(xa)
 
for i in range(N):
    for j in range(N):
        x, y = xa[i,j], ya[i,j]
        if est_ce_conducteur_central(x,y) : 
            V[i,j] = 30
####  Inutile, car V vaut deja 0 partout à l'état initial
#        if est_ce_conducteur_exterieur(x,y) :
#            V[i,j] =0
 
 
def est_ce_conducteur_exterieur(x,y):
   return (((x-xc2)**2/(R2**2))+((y-yc2)**2/(R2**2)))>=1
 
def est_ce_conducteur_central(x,y):
   return (((x-xc1)**2/(R1**2))+((y-yc1)**2/(R1**2)))<=1
 
def est_dans_domaine_de_calcul(x,y):
   return not est_ce_conducteur_central(x,y) and not est_ce_conducteur_exterieur(x,y)
 
Vnew = V.copy()
 
ecart = 1
k=0 
while ecart > 10**-4:
  k+=1
  for i in range(1,N-1): 
    for j in range(1,N-1):
       if est_dans_domaine_de_calcul(xa[i,j],ya[i,j]):
#          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] )             
          Vnew[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1] ) 
  ecart = np.max(np.abs(V - (Vnew))/np.max(V))
  print(k,'   ',ecart)
  V=Vnew.copy()
 
 
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(a,V[int(N/2),:])
 
plt.show()

Il y a d'autres moyens de gagner encore du temps notamment en évitant le "if". Pour ça il faut au préalable calculer les couples d'indices (i,j) surlesquels tu feras les calculs. Je te laisse déjà voir avec le bout de code là.

**Mainot** · 10/12/2017, 20h08

D'accord, j'ai compris !
Merci pour le temps que vous avez pris pour m'apporter toutes ces aides !

**lg_53** · 11/12/2017, 09h34

En vrai je me dois de t'indiquer la suite car là où tu pèche en performance, c'est que tu n'utilise pas numpy comme il est bon de l'utiliser. Numpy est très très rapide, lorsqu'il faut faire des opérations sur des tableaux entier. Mais pour bénéficier de cette rapidité il ne faut pas les parcourir à l'aide de for, il faut utiliser les notations slicées. Ainsi le code ci-dessous est au moins 100 fois plus rapide (ça me prend une demi seconde)! Même si je fais le calcul finalement sur tous le domaine et que j'impose derrière les valeurs dans les conducteurs à chaque itérations, là puissance de numpy est telle que ca reste beaucoup plus valable.

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
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
from time import time
 
 
t0=time() 
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
 
a, b = np.linspace(-1, 19, N), np.linspace(-10, 10, N)
xa, ya = np.meshgrid(a, b) 
V = np.zeros_like(xa)
 
## Liste des indices (i,j) pour lesquels on est dans le conducteur central et exterieur
## Ceci est calculé dès le début et est stocké en mémoire. On ne passe donc pas notre temps à le recalculer à chaque fois
conduc_centr=np.where(((xa-xc1)**2+(ya-yc1)**2)/(R1**2)<=1)
conduc_ext=np.where(((xa-xc2)**2+(ya-yc2)**2)/(R2**2)>=1)
 
V[ conduc_centr ] = 30
 
Vnew = V.copy()
 
ecart = 1
k=0 
while ecart > 10**-4:
  k+=1
  Vnew[1:-1,1:-1] = 0.25*(V[:-2,1:-1] + V[2:,1:-1] + V [1:-1,:-2] + V[1:-1,2:] )
  Vnew[ conduc_centr ] = 30
  Vnew[ conduc_ext ] = 0
  ecart = np.max(np.abs(V - (Vnew))/np.max(V))
  print(k,'   ',ecart)
  V=Vnew.copy()
 
print("Exec time: ",time()-t0," seconds") 
 
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(a,V[int(N/2),:])
 
plt.show()

**Mainot** · 11/12/2017, 21h48

Ah oui, en effet c'est beaucoup plus rapide, je ne connaissais pas numpy (je vais essayer de me documenter) !

Je ne sais pas si ça a à voir mais je ne comprends pas la ligne suivante :

Code :

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

Vnew[1:-1,1:-1] = 0.25*(V[:-2,1:-1] + V[2:,1:-1] + V [1:-1,:-2] + V[1:-1,2:] )

(ce qu'il y a à l'intérieur des crochets).

Merci encore !

**lg_53** · 11/12/2017, 22h45

Ce sont justement les slices !

Exemple, sur un array 1D, ou tu voudrais calculer la somme de 2 termes consécutifs 2 à 2.

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
N=10
a=np.arange(N)
print(a)
 
#### Methode 1, avec for
b=np.zeros(N-1)
for i in range(N) :
     b[i]=a[i]+a[i+1]
print(b)
 
#### Méthode 2, avec des slices
c=a[:-1]+a[1:]
print(c)
 
### Pour mieux comprendre :
print(a[:-1])
print(a[1:])

Le slice permet d'extraire un sous tableau. Il a la forme start:stop:step. step est facultatif et vaut par défaut 1. Si on ne le spécifie pas on écrit donc juste start:stop.
Là tu veut les éléments de ton tableau pour i=0 jusqu'à i=N-1. Donc au lieu de faire un for, tu écris a[0:N-1], que l'on peut raccourcir par a[:-1] (on peut omettre le 0 au début, et si on met un stop négatif, python interprète directement que c'est en nombre d'élement en partant de la fin de la liste).

De même, quand on écrit a[i+1] pour i=0 jusqu'à i=N-1, alors c'est comme d'écrire a[j] pour j=1 jusqu'à j=N, et donc on écrit a[1:N], que l'on peut abréger a[1:].

Essaye d'autre trucs pour mieux comprendre.
Par exemple :

Code :

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

1
2
print(a[1:-3:2])
print(a[2::2])

Et un bon tuto sur numpy t'expliquera tout ca aussi. En voici un petit en anglais, mais avec des examples de base:
https://www.tutorialspoint.com/numpy...nd_slicing.htm

PS : là mon exemple est en 1D. En 2D tu peux choisir de faire des slices sur l'une ou l'autre des dimensions ou bien encore sur les 2, pour obtenir un sous-tableau du tableau initial

**Mainot** · 13/12/2017, 16h10

D'accord, je comprends mieux, merci !

Programme trop long [Python 2.X]

Calcul scientifique Python

Discussions similaires

Partager

Partager