méthode des trapèzes

Version imprimable

Bonsoir,
je possède un algorithme qui calcule l'intégrale simple à l'aide de la méthode des trapèzes, l'algorithme marche sans difficulté, le voici:

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
	SUBROUTINE QTRAP(FUNC1,A9,B9,S9)
	implicit none
	EXTERNAL FUNC1
	DOUBLE PRECISION A9,B9,S9,FUNC1
	DOUBLE PRECISION EPS,OLDS
	Integer J,JMAX
 
	PARAMETER (EPS=1.E-2, JMAX=20)
	IF((B9-A9).LT.1.E-4) THEN
	S9 = 1.E-6
	RETURN
	ENDIF
	IF(B9.EQ.0.0) THEN
	B9 = 0.001
	ENDIF
	OLDS=-1.E30
	DO J=1,JMAX
	CALL TRAPZD(FUNC1,A9,B9,S9,J)
	IF (ABS(S9-OLDS).LT.(EPS*ABS(OLDS))) RETURN
	OLDS=S9
	end do
	RETURN
	END
 
!    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
	SUBROUTINE TRAPZD(FUNC2,A8,B8,S8,N8)
	implicit none
	EXTERNAL FUNC2
	DOUBLE PRECISION A8,B8,S8
	INTEGER N8
	DOUBLE PRECISION SUM,DEL,X,FUNC2
	INTEGER J,IT,TNM
 
 
	!    This routine computes the N'th stage of refinement of an extended
	!    trapezoidal rule. FUNC is input as the name of the function to be
	!    integrated between limits A and B, also input. When called with N=1
	!    the routine returns as S the crudest estimate of the integral.
	!    Subsequent calls with N=2,3,.. will improve the accuracy of S by 
	!    adding 2**N-2 additional interior points. S should not be modified 
	!    BETWEEN SEQUENTIAL CALLS.
 
 
	IF (N8.EQ.1) THEN
	S8=0.5*(B8-A8)*(FUNC2(A8)+FUNC2(B8))
	IT = 1
 
	ELSE
 
	TNM=IT
	DEL=(B8-A8)/TNM
	X=A8+0.5*DEL
	SUM=0.
	DO J=1,IT
			SUM=SUM+FUNC2(X)
			X=X+DEL
	end do
	S8=0.5*(S8+(B8-A8)*SUM/TNM)
	IT=2*IT
	ENDIF
	RETURN
	END

Le problème est dans la variable locale IT, cette variable est déclarée comme un entier mais on ne lui a affecté aucune valeur et pourtant on l'utilise au sein du programme

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
SUBROUTINE TRAPZD(FUNC2,A8,B8,S8,N8)
	implicit none
	EXTERNAL FUNC2
	DOUBLE PRECISION A8,B8,S8
	INTEGER N8
	DOUBLE PRECISION SUM,DEL,X,FUNC2
	INTEGER J,IT,TNM


	!    This routine computes the N'th stage of refinement of an extended
	!    trapezoidal rule. FUNC is input as the name of the function to be
	!    integrated between limits A and B, also input. When called with N=1
	!    the routine returns as S the crudest estimate of the integral.
	!    Subsequent calls with N=2,3,.. will improve the accuracy of S by 
	!    adding 2**N-2 additional interior points. S should not be modified 
	!    BETWEEN SEQUENTIAL CALLS.
    
      
	IF (N8.EQ.1) THEN
	S8=0.5*(B8-A8)*(FUNC2(A8)+FUNC2(B8))
	IT = 1
      
	ELSE
        
	TNM=IT
	DEL=(B8-A8)/TNM
	X=A8+0.5*DEL
	SUM=0.
	DO J=1,IT
			SUM=SUM+FUNC2(X)
			X=X+DEL
	end do
	S8=0.5*(S8+(B8-A8)*SUM/TNM)
	IT=2*IT
	ENDIF
	RETURN
	END

svp aidez moi à comprendre quelles sont les valeurs que prend IT et comment fortran n'affiche aucune erreur malgré le fait qu'on n'a affecté aucune valeur à IT au début
merci de m'aider

14/03/2010, 21h01
FR119492

Salut!
Lorsqu'on écrit une bibliothèque de sous-programmes, il est indispensable de mettre au début de chacun des commentaires indiquant la signification de chaque variable.
Jean-Marc Blanc

bonsoir,
j'ai signalé que l'algorithme calcule l'intégrale simple par la méthode des trapèzes, je n'ai pas d'autres explications, c'est un algorithme préprogrammé mais si vous voulez voila le programme complet:

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
External h0
double precision h0
DOUBLE PRECISION t,f0,x0,a,sigma1,beta1,sigma2,beta2,w,s,x1,x2,w1
COMMON/fonction/ sigma1,beta1,sigma2,beta2,n,i,t
 
CALL QTRAP(h0,0.0,2.0,f0)
print *, f0
end	
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!    LA FONCTION h0
 
DOUBLE PRECISION  FUNCTION h0(a)
DOUBLE PRECISION a,s1,b1,s2,b2,t00,a00
integer n11,i11
real DE
dimension DE(501)
COMMON/fonction/ s1,b1,s2,b2,n11,i11,t00 
h0=a00*(b1/s1)*((a/s1)**(b1-1))*exp(-(a/s1)**b1)
RETURN
END
SUBROUTINE QTRAP(FUNC1,A9,B9,S9)
	implicit none
	EXTERNAL FUNC1
	DOUBLE PRECISION A9,B9,S9,FUNC1
	DOUBLE PRECISION EPS,OLDS
	Integer J,JMAX
 
	PARAMETER (EPS=1.E-2, JMAX=20)
	IF((B9-A9).LT.1.E-4) THEN
	S9 = 1.E-6
	RETURN
	ENDIF
	IF(B9.EQ.0.0) THEN
	B9 = 0.001
	ENDIF
	OLDS=-1.E30
	DO J=1,JMAX
	CALL TRAPZD(FUNC1,A9,B9,S9,J)
	IF (ABS(S9-OLDS).LT.(EPS*ABS(OLDS))) RETURN
	OLDS=S9
	end do
	RETURN
	END
 
!    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
	SUBROUTINE TRAPZD(FUNC2,A8,B8,S8,N8)
	implicit none
	EXTERNAL FUNC2
	DOUBLE PRECISION A8,B8,S8
	INTEGER N8
	DOUBLE PRECISION SUM,DEL,X,FUNC2
	INTEGER J,IT,TNM
 
 
	!    This routine computes the N'th stage of refinement of an extended
	!    trapezoidal rule. FUNC is input as the name of the function to be
	!    integrated between limits A and B, also input. When called with N=1
	!    the routine returns as S the crudest estimate of the integral.
	!    Subsequent calls with N=2,3,.. will improve the accuracy of S by 
	!    adding 2**N-2 additional interior points. S should not be modified 
	!    BETWEEN SEQUENTIAL CALLS.
 
 
	IF (N8.EQ.1) THEN
	S8=0.5*(B8-A8)*(FUNC2(A8)+FUNC2(B8))
	IT = 1
 
	ELSE
 
	TNM=IT
	DEL=(B8-A8)/TNM
	X=A8+0.5*DEL
	SUM=0.
	DO J=1,IT
			SUM=SUM+FUNC2(X)
			X=X+DEL
	end do
	S8=0.5*(S8+(B8-A8)*SUM/TNM)
	IT=2*IT
	ENDIF
	RETURN
	END

15/03/2010, 15h12
Sylvain Bergeron

J'ai regardé rapidement.

IT est le nombre de trapèze. Le calcul est itératif et IT est initialisé à 1 au premier appel. Aux appels suivants, sa valeur est doublée.

Le code se fie au fait que la variable ne perd pas sa valeur d'un appel à l'autre. Ça fonctionne généralement (d'ailleurs, ton programme fonctionne), mais ce n'est pas « standard ». Ce comportement est géré par l'attribut SAVE dans le langage, mais cet attribut est souvent implicite. Tu devrais le mettre explicitement pour que le programme soit plus facile à comprendre (et que le résultat ne dépendent pas de certaines switches de compilation).