Numpy et fft

**lafrite26** · 27/06/2010, 20h43

Bonjour,

Depuis une semaine je me casse la tête sur le module fft de numpy. Et je ne comprend pas comment il marche!
Je suis rendu compte que fft.fft et fft.ifft n'était pas l'inverse l'une de l'autre!
Si je lance ce script (pas très propre, j'ai juste extrait quelques lignes d'un plus gros code pour tester fft):

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
60
61
62
 
from math import *
import numpy as np
 
# Boundaries
xMin = 0
xMax = 2*pi
 
# Number of points
N = 256
h = (xMax - xMin)/N
 
tMax = 5.
numberLaps = 20
dt = tMax/numberLaps
 
x = h*np.arange(N) - (xMax - xMax)/2
 
t = 0
 
v=[]
for xI in x:
	v += [exp(-10*(xI-pi)**2)]
 
v = np.array(v)
vOld = v
 
diffMax = []
diffMin = []
 
for i in np.arange(1,numberLaps):
	t = t + dt
	print '############ t = ' , t
 
	v_hat = np.fft.fft(v)
 
	vBis = np.fft.ifft(v_hat) 
 
	print np.allclose(v.astype(complex), np.fft.ifft(np.fft.fft(v)))
	print np.allclose(v.astype(complex), vBis)
 
 
	diffMax += [(v - vBis.real).max()]
	print 'diff max' , (v - vBis.real).max()
	diffMin += [(v - vBis.real).min()]
	print 'diff min' , (v - vBis.real).min()
 
	list_k = np.array(range(v_hat.shape[0]/2+1) + range (-v_hat.shape[0]/2+1 , 0))
	vTT_hat = - (list_k)**2 * v_hat
 
	D2v =  np.fft.ifft(vTT_hat).real
 
	vNew = dt**2 * D2v + 2*v - vOld
 
	vOld = v
	v = vNew
 
	print ''
 
 
print 'diff max' , max(diffMax)
print 'diff min' , max(diffMin)

Au début la difference entre la "vrai" fonction et celle transformé puis retransformé est correcte (10^-16). Mais au petit à petit cette difference augmente pour atteindre 10^22).
Mais chose encore plus étrange, allclose me dit qu'elles sont effectivement proches!
Je vous colle ce qu'il le script renvoie:

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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
 
############ t =  0.25
True
True
diff max 1.6555156296e-16
diff min -1.38777878078e-16
 
############ t =  0.5
True
True
diff max 1.52655665961e-16
diff min -1.11022302463e-16
 
############ t =  0.75
True
True
diff max 2.74086309204e-16
diff min -2.42912858918e-16
 
############ t =  1.0
True
True
diff max 1.86243314806e-15
diff min -1.7763568394e-15
 
############ t =  1.25
True
True
diff max 1.73702300102e-14
diff min -1.57884772964e-14
 
############ t =  1.5
True
True
diff max 2.09339490187e-13
diff min -2.30038210702e-13
 
############ t =  1.75
True
True
diff max 2.04636307899e-12
diff min -2.33058017329e-12
 
############ t =  2.0
True
True
diff max 5.45696821064e-11
diff min -7.27595761418e-11
 
############ t =  2.25
True
True
diff max 4.47034835815e-08
diff min -3.72529029846e-08
 
############ t =  2.5
True
True
diff max 3.43322753906e-05
diff min -4.57763671875e-05
 
############ t =  2.75
True
True
diff max 0.03515625
diff min -0.03515625
 
############ t =  3.0
True
True
diff max 34.0
diff min -40.0
 
############ t =  3.25
True
True
diff max 36864.0
diff min -36864.0
 
############ t =  3.5
True
True
diff max 37748736.0
diff min -37748736.0
 
############ t =  3.75
True
True
diff max 38654705664.0
diff min -34359738368.0
 
############ t =  4.0
True
True
diff max 3.51843720888e+13
diff min -3.95824185999e+13
 
############ t =  4.25
True
True
diff max 3.6028797019e+16
diff min -3.6028797019e+16
 
############ t =  4.5
True
True
diff max 3.45876451382e+19
diff min -3.68934881474e+19
 
############ t =  4.75
True
True
diff max 3.7778931863e+22
diff min -3.7778931863e+22
 
diff max 3.7778931863e+22
diff min -1.11022302463e-16

Si quelqu'un peut m'expliquer pourquoi est ce que j'ai ce comportement? Jai certainement louper quelque chose dans la manipulation de fft.
Merci d'avance

**Matthieu Brucher** · 27/06/2010, 23h48

Si, ifft(fft()) retourne aux erreurs numériques près, les mêmes valeurs.
Plusieurs points :
- tu utilises très insiffisamment numpy, l'initialisation de v est réalisable en une seule opération, bien meilleur niveau efficacité
- ton erreur doit provenir des calculs que tu effectues, mais je n'ai pas le temps de déchiffrer ce que tu fais pour savoir quel est l'ago sous-jacent...

**lafrite26** · 28/06/2010, 01h53

Ok tant mieux ça me rassure que fft(ifft) soit bien l'identité!
Mais quand même, comment ça se fait que j'arrive avec des écarts de l'ordre de 10^22?
Bien que je le fasse l'un près l'autre de la manière suivante,

Code :

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

1
2
3
4
 
	v_hat = np.fft.fft(v)
	vBis = np.fft.ifft(v_hat) 
print (v - vBis.real).max()

L'écart grandit rapidement au bout de quelques iterations (20 dans ce programme!).

Voici une version plus épuré avec des commentaires de ce que j'aimerai faire tourner correctement. Je fais peut être un grosse bétise sans m'en rendre compte.

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
 
# -*- coding: utf-8 -*-
from math import *
import numpy as np
 
# interval
xMin = 0
xMax = 2*pi
 
# Nombre de points
N = 16
 
# Pas
h = (xMax - xMin)/N
 
# Temps
tMax = 5.
numberLaps = 200
dt = tMax/numberLaps
 
# Points
x = h*np.arange(N)
 
# Initialisation de v et vOld
v=np.exp(-10*(xI-pi)**2)
vOld = v
 
# Stockage
dataV = [v]
 
# Recurence sur le temps
for i in np.arange(1,numberLaps):
 
	# On passe dans le monde Fourier
	v_hat = np.fft.fft(v)
 
	# On derive deux fois v dans le monde Fourier
	list_k = np.array(range(v_hat.shape[0]/2+1) + range (-v_hat.shape[0]/2+1 , 0))
	vTT_hat = - (list_k)**2 * v_hat
 
	# On revient dans le monde réel
	D2v =  np.fft.ifft(vTT_hat).real
 
	# On calcul le nouveau v (l'équation differencielle est u_tt - u_xx = 0)
	vNew = dt**2 * D2v + 2*v - vOld
 
	# On prépare pour la prochaine étape
	vOld = v
	v = vNew
 
	# On stock le nouveau v
	dataV += [v]

Concernant l'initialisation de v, elle était écrite comme ça car avant j'avais testé ce programme avec une fonction définit par morceau et que je ne sais pas faire autrement pour ce type de fonction. Mais la version de ce post devrait plus te convenir?

Merci pour ton aide

**Matthieu Brucher** · 28/06/2010, 09h19

Si au début, c'est bon, c'est que tu as un problème avec ton calcul qui n'est pas bon.

**lafrite26** · 28/06/2010, 15h24

Je ne comprend pas où je peux me planter quand je teste si la fonction et la fonction transformée puis retransformée sont les même à travers ce bloc d'instruction?

Code :

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

1
2
3
4
 
	v_hat = np.fft.fft(v)
	vBis = np.fft.ifft(v_hat) 
print (v - vBis.real).max()

Normalement je devrais voir afficher des valeurs proches de 0 non?

**Matthieu Brucher** · 28/06/2010, 16h16

De rien... Quelle est ta version de Numpy ?

**lafrite26** · 28/06/2010, 21h29

Heu non en fait, jai supprimé le post où je disais que le problème venait d'autre peut être de la version de numpy installé sur mon ordi.
Donc le problème reste entier. Je ne comprend pas pourquoi j'ai des valeurs très grande quand j'execute ça dans mon programme?

Code :

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

1
2
3
4
 
	v_hat = np.fft.fft(v)
	vBis = np.fft.ifft(v_hat) 
print (v - vBis.real).max()

Numpy et fft

Calcul scientifique Python

Discussions similaires

Partager

Partager