Discussion :

[lisoux]

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