Numpy et fft

Version imprimable

Numpy et fft

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:

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:

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

27/06/2010, 22h48
Matthieu Brucher

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...

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:

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:

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

28/06/2010, 08h19
Matthieu Brucher

Si au début, c'est bon, c'est que tu as un problème avec ton calcul qui n'est pas bon.
28/06/2010, 14h24
lafrite26
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:

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?
28/06/2010, 15h16
Matthieu Brucher

De rien... Quelle est ta version de Numpy ?
28/06/2010, 20h29
lafrite26
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:

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