intégration sur fortran

Bonsoir
J’ai écrit un programme sur fortran qui calcule l’intégrale d’une fonction h4 à l’aide de la fonction QTRAP, le problème est que ce programme ne génère aucune erreur mais donne des résultats aberrants :1.00000E-6 pour toutes les valeurs de i et même si je change ma fonction à intégrer les mêmes résultats apparaissent.
Voici le 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
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
program verif h4
	EXTERNAL h4
 
      DOUBLE PRECISION t,sigma1,beta1,w
	integer n,i
 
	COMMON/fonction/ sigma1,beta1,n,i,t
 
      real f
	dimension f(5001)
 
	sigma1=1.0
	beta1=2.0
      t=5.0
	n=t*500.0
 
	DO i=1,n,1
	CALL QTRAP(h4,0.0,2.0,w)
      f(i)=w
	print*,f(i)
	ENDDO
	end
 
!    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!    LA FONCTION à intégrer
!    ALPHA: PARAMETRE DE FORME
!    PLMD: PARAMETRE D'ECHELLE 	
		DOUBLE PRECISION  FUNCTION h4(a)
	 	DOUBLE PRECISION a,s1,b1,t00
	    integer i11,n11
	    COMMON/fonction/ s1,b1,n11,i11,t00
 
	h4=1-(1-exp(-((i11-a)/s1)**b1))
 
 		RETURN
		END
 
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
		SUBROUTINE QTRAP(FUNC1,A9,B9,S9)
		EXTERNAL FUNC1
		DOUBLE PRECISION A9,B9,S9
		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)
		EXTERNAL FUNC2
		DOUBLE PRECISION A8,B8,S8
		INTEGER N8
		DOUBLE PRECISION SUM,DEL,X
		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 car je n’ai pas pu trouver l’origine de ce problème car je vois que mes variables sont bien déclarés.
merci

Salut!

Citation:

Svp aidez moi car je n’ai pas pu trouver l’origine de ce problème

Ajoute

Code:

Implicit None

au début de ton programme principal et de chacun de tes sous-programmes (c'est une règle à appliquer systématiquement pour que les erreurs apparaissent déjà à la compilation).

Ensuite, il y a une méthode générale à appliquer quand on se trouve dans ton cas: saupoudre des Write un peu partout dans ton code et regarde où ça commence à devenir faux.
Jean-Marc Blanc

bonjour,
merci bcp pour la réponse, j'ai mis la implicit none devant chaque programme et sous programme et l'erreur suivante persiste à apparaitre:
Code:

1 2 Error: This name does not have a type, and must have an explicit type. [FUNC2] S8=0.5*(B8-A8)*(FUNC2(A8)+FUNC2(B8))
sachant que S8 est bel et bien déclarée comme un réel et FUNC2 est déclarée comme un external

svp aidez moi
merci

Salut!
Comment veux-tu que ton ordinateur sache que FUNC2 est de type double précision?
Jean-Marc Blanc

bonjour,
merci beaucoup, votre aide m'a été précieuse.