Boucles imbriquées, intégrale double, tableaux numpy

Salut !


J'ai le code suivant, dans lequel je n'arrive pas à remplacer la double boucle imbriquée par une opération sur des tableaux numpy. La double boucle est beaucoup trop lente (je dois appeler la fonction des centaines de fois et les tableaux font 1000x800 :

Code:

---

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

def GetFlux(self, time):   bx = self.GetField("bx", time) * self.wpewce by = self.GetField("by", time) * self.wpewce bz = self.GetField("bz", time) * self.wpewce       flux = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN') flux2 = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN')   dx = self.dl[0] dz = self.dl[1]   nx = self.ncells[0] nz = self.ncells[1]   j = 0   for i in np.arange(1, nx): flux2[i,0] = flux2[i-1,0] + bz[i-1,0]*dx   flux[1:,0] = flux[0,0] + np.cumsum(bz[:-1,0]*dx)       for j in np.arange(1,nz): flux2[0,j] = flux2[0,j-1] - bx[0,j-1]*dz   flux[0,1:] = flux[0,0] - np.cumsum(bx[0,:-1]*dz)     for i in np.arange(1,nx): for j in np.arange(1,nz): flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)     # flux = ???     return flux2

---

(voir la me chose ici : http://pastebin.com/0ZEr6hKL)


Le but est de calculer l'intégrale double donnant F a partir de F=dB_z/dx et F=-dB_x/dz.


'flux2' est calculé à partir de simples boucles, comme on pourrait le faire en C, et contient le bon résultat.

'flux' devrait contenir la meme chose, mais fait a partir d'opérations sur tableaux avec un fancy indexing. Seulement je ne parviens pas a faire la double boucle.


Attention, c'est une intégrale, donc dans les boucles;, c'est bien flux2[i,j] qui dépend de flux2[i-1,j] et flux2[i,j-1], i.e. le tableau dépend de lui même à litération d'avant... c'est ce qui me pose problème pour la boucle 2D.

Dans les boucles simples plus haut, j'ai trouvé cumsum() pour bien faire le travail, et les cases de tableau concernées sont égales à celles obtenues dans les boucles simples.

Je me dis qu'il doit exister une méthode basée sur cumsum pour la boucle imbriquée ?


Merci d'avance bcp pour votre aide :)

Bonjour,

Avez-vous pensez à un produit de convolution ( convolve) avec dans ce cas un masque 3x3 pout flux2 3x1 et 1x3 pour bz et bx.

Bonjour.

Je n'ai pas compris comment Heimdall a obtenu l'équation sur flux2. Alors j'ai déterminé le flux en intégrant l'équation de Maxwell-Ampère.
En discrétisant, on montre que le flux est proportionnel à :

Code:

---

dz * np.cumsum( bz, axis=0) - dx * np.cumsum( bx, axis=1)

---

Cette expression est bien évidemment différente de celle donnée pour calculer flux2.

Pour les comparer, je fais l'hypothèse que le champ magnétique est suivant X (Bz=0).
Je modélise le champ par une gaussienne 2D, de révolution autour de l'axe Y.
La gaussienne est placée au centre du domaine. Soient larg l'écart type, et maxi la valeur maximale des abscisses. Je note r le rapport larg sur maxi.

Pour r=3, la gaussienne est "étendue" par rapport au domaine, le champ tend à être uniforme (Fig. ChampPresqueUniforme).
Lorsque le champ est rigoureusement uniforme, les résultats sont les mêmes : Bx a une valeur constante dans les plans perpendiculaires à l'axe X.
Un léger décalage apparait à r=3 avec la formule de Heimdall vers les Z>0. Avec ma formule, les courbes d'iso-valeurs sont de révolution autour de l'axe X.

Pour r=1/3, la gaussienne est peu "étendue" sur le domaine, avec un comportement asymptotique de champ nul. (Fig. ChampPresqueNul).
Avec ma méthode les courbes d'iso-valeurs restent de révolution autour de X, tandis qu'avec celle de Heimdall, elles sont clairement orientées suivant un axe faisant un angle de 45° avec X.

Qui a raison ?
Quelqu'un connaitrait-il un cas théorique permettant de trancher ?