[débutante] problème algorithme intégrale simpson et séries Fourier

**lisoux** · 26/12/2014, 17h47

Salut à tous,

J'ai un problème dans ma sous-routine de calcul d'intégrale par la méthode de Simpson. Lorsque j'augmente le nombre de subdivisions, le résultat devient de moins en moins précis alors que ça devrait être le contraire. J'ai demandé l'aide de mon prof qui ne m'a jamais répondu.

Voici mon code:

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
subroutine simpson (a,b,n,l)
!calcul de l'intégrale par la méthode de simspon
!pour simpson il faut que n soit pair!
 
use kinds
 
implicit none
 
!a= borne inférieure
!b = borne supérieure
!n= nombre de divisions
!l= résultat de l'intégration par simpson
 
real(kind=r),intent(in)::a,b
integer,intent(in)::n
real(kind=r),intent(out)::l
 
!h= pas d'intégration
!x0=a (borne inférieure)
!xn=b (borne supérieure)
!u est la valeur de la première somme * le coefficient devant(celle allant de 1 à n/2-1)
!v est la valeur de la deuxième somme * le coefficient devant (celle allant de 1 à n/2)
 
real(kind=r)::h,x0,xn,u,v
 
!i,j= indices boucles
!q,s= indices des bornes supérieures des sommes
 
integer::i,j,q,s
!x= vecteur contenant les x(2i)
!z=vecteur contenant les valeurs de sin(x(2i))=f(x(2i))
!y= vecteur contenant les x(2i-1)
!d=vecteur contenant les valeurs de sin(x(2i-1))=f(x(2i-1))
 
real(kind=r),allocatable,dimension(:)::x
real(kind=r),allocatable,dimension(:)::z
real(kind=r),allocatable,dimension(:)::d
real(kind=r),allocatable,dimension(:)::y
 
!q,s bornes sup sommes
 
q=int(n/2.)	
s=int(n/2.-1.)
 
!j'alloue les tailles aux vecteurs
 
allocate(y(1:q))	
allocate(z(1:s))
allocate(x(1:s))
allocate(d(1:q))
 
!je définis x0 et xn
 
x0=a
xn=b
 
!je calcule le pas d'intégration
 
h=real((b-a)/real(n))
 
	!boucle pour la première somme
	 do i=1,s
 
 
		x(i)=a+2.*real(i)*h
 
		z(i)=sin(x(i))
 
	end do	
 
	u=2.*sum(z)
 
	!boucle pour la 2ème somme
	do j=1,q
 
 		y(j)=a+(2.*real(i)-1.)*h
 
 		d(j)=sin(y(j))
 
 	end do
 
 	v=4.*sum(d)
 !valeurs de l'intégrale par simpson
 
 l=(h/3.)*(sin(x0)+v+u+sin(xn))
 
 
 deallocate(y)	
 deallocate(z)
 deallocate(x)
 deallocate(d)
 
end subroutine simpson

J'ai un second problème. Dans un autre programme( le calcul des coefficients des séries de Fourier), j'ai une fonction impaire donc les coefficients a0 et an doivent valoir 0. Mon a0 vaut 0 mais pas le an. Et je pense que le bn n'est pas bon non plus.

voici le programme principal:

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
!*******************************************************************************
! program fourier
! programme calculant la série de fourier d'une fonction et affichant sur un 
! graphique la fonction et la série
! 
!*******************************************************************************
!les a valent 0 car fonction impaire
program fourier
 
use kinds
use interface1
use interface2
use constante
implicit none
 
real(kind=r)::x,L
real(kind=r)::a,b,a0,an,bn
integer::n
real(kind=r)::m
integer::harm,i,j
!real(kind=r),allocatable,dimension(:)::ser
!real(kind=r),allocatable,dimension(:)::y
 
 
 
 
write(*,*)"entrez la periode L"
read(*,*)L
 
 open(12,file='parameter.dat')
 
     write(12,*) L
 
 close (12)
 
 !allocate(ser(0:int(L)))
 !allocate(y(0:int(L)))
a=0.
n=100.
!calcul de a0
	call trapeze1(a,L,n,m)
	a0=(1./L)*m
 
	write(*,'(a,f15.10)')"a0=",a0
 
	call trapeze2(a,L,n,m)
	!write(*,*)m
	an=(2./L)*m
 
	write(*,'(a,f15.10)')"an=",an
 
	call trapeze3(a,L,n,m)
	bn=(2./L)*m
 
	write(*,'(a,f15.10)')"bn=",bn
 
write(*,*)"entrez le nombre d'harmonique souhaite"
read(*,*)harm
 
!do i=0,int(L)
 
	!y(i)=real(i)
	!ser(i)=serie(y(i),harm,bn)
 
 
!end do
!write(*,*)y,ser
 
 
call PGBEGIN (0, 'pgplot.ps/cps', 1, 1) !ouvrir un dispositif graphique
 ! Unit,type de device (ps pour postscript), nombre de subdivisions
 call PGSLW(3)                           ! choisir la largeur des traits
 call PGSCH(1.5)                         ! choisir la hauteur des caracteres
 call PGPAGE                             ! accès à une nouvelle subdivision
 call PGVSTD                             ! definir la viewport "standard"
 call PGBOX('BCNST',0.0,0,'BCNST',0.0,0) ! choisir le style des axes
 call PGSWIN(-0.5,L+1,-1.1,1.1)        ! choisir le système de coordonnées
 call PGLABEL('x', 'SF', 'série de fourier') ! légendes et titre
 call PGSCI (4)                          ! choisir la couleur (bleu)
!call PGPOINT(n, y, ser, 4)              ! tracer des points
 !call PGSCI (2)                         ! choisir la couleur (rouge)
 !call PGSLS (1)                          ! choisir le style des lignes (solide)
 !call PGLINE(n, y, ser)               ! relie par des segments de droite 
 !call PGSCI(6)
 call PGFUNX (fonction,n,0,L,1)
 !call PGSCI (2)
 !call PGFUNX(serie,n,0,L,1)
 
 
 
 call PGEND				 ! terminer le graphique
 
!deallocate(ser)	
!deallocate(y)
 
end program

Voici la sous-routine qui calcule les an

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
subroutine trapeze2 (a,b,n,m)
 
!calcul de l'intégrale par la méthode des trapèzes
use kinds
use interface1
use constante
implicit none
 
!a= borne inférieure
!b = borne supérieure
!n= nombre de divisions
!m= résultat de l'intégration par la méthode des trapèzes
 
real(kind=r),intent(in)::a,b
integer,intent(in)::n
real(kind=r),intent(out)::m
 
!h= pas d'intégration
!y est la première partie de l'expression de l'intégrale
 
real(kind=r)::h,y,L
 
 
!x= vecteur contenant les valeurs des x(i)
!z= vecteur contenant les valeurs de h*f(x(i))
 
real(kind=r),allocatable,dimension(:)::x
real(kind=r),allocatable,dimension(:)::z
 
!indice boucle
 
integer::i
 
!j'alloue une taille aux vecteurs
 
allocate(x(2:n-1))	
allocate(z(2:n-1))
 
!je calcule le pas d'intégration
 
h=real((b-a)/real(n))
 
open(12,file='parameter.dat')
read(12,*)L
 
 
!je calcule le premier terme de l'expression de l'intégrale
 
y=((fonction(a)*cos((2.*pi*real(n)*a)/L)+fonction(L)*cos((2.*pi*real(n)*L)/L))/2.)*h
 
	!boucle calculant (xi) et f(x(i))
 
 
 
 		do i=2,n-1
 
 
				x(i)=a+real(i)*h
 
 
			z(i)=h*(fonction(x(i))*cos((2.*pi*real(n)*x(i))/L))
 
	!valeurs de l'intégrale par la méthode des trapèzes
 
	m=sum(z)+y
 
end do	
 deallocate(x)
 deallocate(z)
 
close(12)	
end subroutine trapeze2

Voici la sous-routine qui calcule les bn (très peu de changements par rapport à la sous-routine précédente):

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
subroutine trapeze3 (a,b,n,m)
 
!calcul de l'intégrale par la méthode des trapèzes
use kinds
use interface1
use constante
implicit none
 
!a= borne inférieure
!b = borne supérieure
!n= nombre de divisions
!m= résultat de l'intégration par la méthode des trapèzes
 
real(kind=r),intent(in)::a,b
integer,intent(in)::n
real(kind=r),intent(out)::m
 
!h= pas d'intégration
!y est la première partie de l'expression de l'intégrale
!x0=a=borne inférieure d'intégration
!xn=b=borne supérieure d'intégration
 
real(kind=r)::h,y,x0,xn,L
 
 
!x= vecteur contenant les valeurs des x(i)
!z= vecteur contenant les valeurs de h*f(x(i))
 
real(kind=r),allocatable,dimension(:)::x
real(kind=r),allocatable,dimension(:)::z
 
!indice boucle
 
integer::i
 
!j'alloue une taille aux vecteurs
 
allocate(x(2:n-1))	
allocate(z(2:n-1))
 
!je calcule le pas d'intégration
 
h=real((b-a)/real(n))
 
open(12,file='parameter.dat')
read(12,*)L
 
!je calcule le premier terme de l'expresison de l'intégrale
 
y=((fonction(a)*sin((2.*pi*real(n)*a)/L)+fonction(L)*sin((2.*pi*real(n)*L)/L))/2.)*h
 
	!boucle calculant (xi) et f(x(i))
 
 	do i=2,n-1
 
 
			x(i)=a+real(i)*h
 
 
		z(i)=h*fonction(x(i))*sin((2.*pi*real(n)*x(i))/L)
	end do	
	!valeurs de l'intégrale par la méthode des trapèzes
 
	m=sum(z)+y
 deallocate(x)
 deallocate(z)
close(12)	
end subroutine trapeze3

Voici la fonction fournie par le prof: dans celle-ci j'ai également un problème: lorsque j'essaie de lire dans le ficher la variable L, ca me dit "end of file" or j'ai vérifié partout si j'avais bien refermé le fichier une fois que je l'avais lu (sinon ca aurait pu expliquer l'erreur)

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
function fonction (x)
 
use kinds
 
implicit none
 
real(kind=r),intent(in):: x
 
real(kind=r)::fonction,L 
 
open(12,file='parameter.dat')
 
read(12,*)L
 
 
if (0<x .and. x<(L/2.)) then 
 
	fonction=+1.
 
elseif (L/2.<x .and. x<L) then
 
	fonction=-1.
 
else if (x==0. .or. x==L/2.) then
	fonction=0.
end if
 
 
close(12)
end function fonction

et voici la série de fourier obtenue en utilisant le coefficient bn:

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
function serie(x,harm,bn)
 
use kinds
use constante
implicit none
 
real(kind=r),intent(in)::x,bn
integer,intent(in)::harm
real(kind=r)::serie,L
real(kind=r),allocatable,dimension(:)::S
integer::i
 
 
open(12,file='parameter.dat')
read(12,*)L
 
allocate(S(1:harm))
 
		do i=1,harm
 
		!y(i)=real(i) 
 
 
			S(i)=bn*sin((2.*pi*real(i)*x)/L)
 
		end do
		serie=sum(S)
 
close(12)	
end function serie

Merci beaucoup d'avance et bonne soirée

**Aotarie** · 28/12/2014, 22h55

Bonsoir,
Sans savoir si le problème va être résolut, une erreur d'indice est présente ligne 76 dans la sub de simpson : un "j" irait mieux qu'un "i".(en fait je crois que c'est là d’où vient la faute

)
Par ailleurs, tu te compliques mais alors beaucoup la vie (les tableaux, moi aussi j'aime ça, mais faut se limiter). J'ai repris ton algo (ça va me servir pour plus tard) en allant de façon plus directe :

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
 
        subroutine integrale(a,b,n,l)
		real,intent(in)::a,b
		real,intent(out)::l
		real::h,x,y
		integer,intent(in)::n
 
		x=0
		y=0
		h=real((b-a)/real(n))
 
		do i=1,n/2-1
			x=sin(a+i*2*h)+x
		end do
 
		do j=1,n/2
			y=sin(a+(2*j-1)*h)+y
		end do
 
		l=(h/3)*(sin(a) + 2*x + 4*y + sin(b))
 
	end subroutine

Si jamais tu fais des testes, faire gaffe aux unités (rad, deg, grad) sinon on cherche longtemps "pk qui veut pas marcher cet algo

! ! ! "

Cordialement

**lisoux** · 28/12/2014, 23h23

Bonsoir,

Ce programme marche enfin à merveille, je vous remercie énormément

Il ne reste plus qu'à trouver une personne aussi gentille que vous pour trouver l'erreur dans les séries de Fourier car ce programme me désespère.

Bonne soirée

Bien à vous

**Aotarie** · 29/12/2014, 14h23

Bonjour,
Je n'ai pas encore trouvé pourquoi cela ne marchait pas mais déjà pour y voir plus clair, voici 2-3 trucs :

Avec ça tu définis tes variables en globales, donc tu peux laisser les a,b,L tout au long du programme.

Code :

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

1
2
3
4
5
 
Module variable
	implicit none
	real::L,a,b
end module variable

Ensuite, faut bien les définir quelque part donc :

Code :

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

1
2
3
4
5
6
7
8
9
10
 
subroutine chargement_variable()
	use variable
	implicit none
	open(11,file='L.dat')
	read(11,*)L
	close(11)
        a=0
	b=L
end subroutine

Après, ça serait cool si la fonction f était définit une bonne fois pour toute :

Code :

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

1
2
3
4
5
6
7
 
real function f(x)
	use variable		
	if(x>0 .and. x<L/2)f=1
	if(x>L/2 .and. x<L)f=-1	
	if(x==0 .or. x==L/2)f=0
end function f

Pareil, si tu ne veux pas te trimbaler les cos et sin :

Code :

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

1
2
3
4
5
 
real function g(x,n)
	use variable
	g=(cos(2*3.14/L*n*x))
end function g

et :

Code :

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

1
2
3
4
5
 
real function t(x,n)
	use variable
	t=(sin(2*3.14/L*n*x))
end function t

Ensuite, une sub d'intégration, c'est beau (elle est belle ma sub, non ?

) mais une fonction me parait plus facile à utiliser ici, du coup :

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
 
real function integration(borne_inf,borne_sup,qualite)	
	real::h,tmp1,tmp2,x
	integer,intent(in)::qualite
	tmp1=0
	tmp2=0
	h=real((borne_sup-borne_inf)/real(qualite))
	do i=1,qualite/2-1
		x=(borne_inf+i*2*h)
		tmp1=f(x)+tmp1
	end do
	do j=1,qualite/2
		x=(borne_inf+(2*j-1)*h)
		tmp2=f(x)+tmp2
	end do
	integration=(h/3)*(f(borne_inf) + 2*tmp1 + 4*tmp2 + f(borne_sup))
end function integration

Avec ça, y a moyen de le finir plus rapidement ce TD.

Cordialement

**lisoux** · 29/12/2014, 14h46

Salut

Effectivement avec tes conseils, le programme se simplifie pas mal...

Merci

**Aotarie** · 29/12/2014, 14h57

Pourras tu afficher ton code à la fin, stp ? (s'il finit par marcher)

**lisoux** · 29/12/2014, 15h27

oui bien sûr

**lisoux** · 30/12/2014, 09h42

Salut!

J'ai une bonne et une moins bonne nouvelle.

La bonne c'est que j'ai remplacé les sous-routines qui faisaient l'intégrale par la méthode des trapèzes par celles de simpson et maintenant j'obtiens bien que a0 et an valent 0 (et donc j'en déduis que le bn est juste).

La moins bonne c'est que dès que j'introduis un graphique pour dessiner la fonction fournie et la série de Fourier qui en résulte, j'obtiens un message d'erreur et je ne vois pas le rapport entre cette erreur et le graphique.

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
call PGBEGIN (0, 'pgplot.ps/cps', 1, 1) !ouvrir un dispositif graphique
 ! Unit,type de device (ps pour postscript), nombre de subdivisions
 call PGSLW(3)                           ! choisir la largeur des traits
 call PGSCH(1.5)                         ! choisir la hauteur des caracteres
 call PGPAGE                             ! accès à une nouvelle subdivision
 call PGVSTD                             ! definir la viewport "standard"
 call PGBOX('BCNST',0.0,0,'BCNST',0.0,0) ! choisir le style des axes
 call PGSWIN(-0.5,L+1,-1.1,1.1)        ! choisir le système de coordonnées
 call PGLABEL('x', 'SF', 'série de fourier') ! légendes et titre
 call PGSCI(6)
 call PGFUNX (fonction,n,0,L,1)
 call PGSCI (2)
 call PGFUNX(serie,n,0,L,1)
 call PGEND				 ! terminer le graphique

Nom : erreur.jpg
Affichages : 980
Taille : 254,3 Ko

**Aotarie** · 30/12/2014, 13h42

Je ne sais pas si cela va changer quelque chose, mais si tu mets tous les modules, fonctions et sub dans le même programme (exemple.f90), le rapport indique t-il la même chose ?

Après, je tenterais une autre approche : j'écrirai mes toutes mes valeurs dans un fichier texte , puis je tracerai à partir de celui-ci.
Un truc du style :
x val_fonction val_fourier (à ne pas mettre dans le fichier .txt)
1 . . . . . .
2
...

Ça a l'avantage de voir les valeurs sorties et de repérer plus rapidement des erreurs. (bon c'est faisable comme il n'y a pas bcp de données)

**lisoux** · 03/01/2015, 16h46

Salut

J'ai fait 3 fichiers: -le programme principal - les sous-routines et fonctions -les modules. Du coup j'ai pu exécuter mais la graphique n'est pas juste. La série de Fourier donne une droite.
Malheureusement, je n'ai plus le temps de m'en occuper

Merci encore

[débutante] problème algorithme intégrale simpson et séries Fourier

Fortran

Discussions similaires

Partager

Partager