double intégrale fortran

**suzanne1307** · 09/07/2009, 14h32

Bonjour,
j'ai utilisé la fonction DGQINT45 pour le calcul du double intégrale d'une fonction exp(y/x), voici le code de la fonction:

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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
 
C***********************************************************************
C                                                                      *
C                  DOUBLE INTEGRAL ALGORITHM 4.5                       *
C                                                                      *
C***********************************************************************
C
C
C
C     TO APPROXIMATE I = DOUBLE INTEGRAL (( F(X,Y) DY DX )) WITH LIMITS
C         OF INTEGRATION FROM A TO B FOR X AND FROM C(X) TO D(X) FOR Y:
C
C     INPUT:   ENDPOINTS A,B; POSITIVE INTEGERS M, N.  (ASSUME THAT THE
C              ROOTS R(I,J) AND COEFFICIENTS C(I,J) ARE AVAILABLE FOR
C              I EQUALS M AND N AND FOR 1<=J<=I.)
C
C     OUTPUT:  APPROXIMATION J TO I.
C
      program double integrale
      implicit none
	double precision XZ,YZ,Y,A,B,H1,H2,AJ,X,C1,D1,Q,JX,K1,K2
	integer M,N,J,I
 
	real R,CO
      DIMENSION R(5,5),CO(5,5)
      real C,D,F
      DIMENSION C(30),D(30),F(5,5)
      CHARACTER*1 AA
      LOGICAL OK
	A=0.1
	B=0.5
	M=5.0
	N=5.0
C     CHANGE FUNCTIONS F, C, D FOR A NEW PROBLEM
C     LIMITS OF INTEGRATION
C     C IS THE LOWER LIMIT OF Y
      C(XZ)=XZ**3
 
C     D IS THE UPPER LIMIT OF Y
      D(XZ)=XZ**2
 
C     DEFINE INTEGRAND FUNCTION F(X,Y)
      !F(XZ,YZ)=X*Y*Z
	F(XZ,YZ)=EXP(YZ/XZ)
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is Gaussian Quadrature for double integrals.'
      WRITE(6,*) 'Have the functions F, C, and D been created?'
      WRITE(6,*) 'Enter Y or N '
      WRITE(6,*) ' '
      READ(5,*)  AA
      IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN
         OK = .FALSE.
10       IF (OK) GOTO 11
            WRITE(6,*) 'Input lower limit of integration and'
            WRITE(6,*) 'upper limit of integration separated'
            WRITE(6,*) 'by a blank.'
            WRITE(6,*) ' '
            READ(5,*) A, B
            IF (A.GE.B) THEN
               WRITE(6,*) 'lower limit must be less than upper limit'
               WRITE(6,*) ' '
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 10
11       OK = .FALSE.
14       IF (OK) GOTO 15
            WRITE(6,*) 'Input two positive integers M, N.'
            WRITE(6,*) 'Both must be less than or equal to 5'
            WRITE(6,*) 'for this implementation.'
            WRITE(6,*) 'Gaussian Quadrature uses M for the outer'
            WRITE(6,*) 'integral and N for the inner'
            WRITE(6,*) 'integral - separate with blank.'
            WRITE(6,*) ' '
            READ(5,*) M,N
            IF(M.LT.6.AND.N.LT.6.AND.N.GT.0.AND.M.GT.0) THEN
              OK=.TRUE.
            ELSE
              WRITE(6,*) 'Both must be positive integers less than'
              WRITE(6,*) 'or equal to 5 '
              WRITE(6,*)
            ENDIF
         GOTO 14
15       CONTINUE
      ELSE
         WRITE(6,*) 'The program will end so that the functions'
         WRITE(6,*) 'F, C and D can be created '
         OK = .FALSE.
      ENDIF
      IF (.NOT.OK) GOTO 400
      R(2,1) = 0.5773502692
      R(2,2) = -R(2,1)
      CO(2,1) = 1.0
      CO(2,2) = 1.0
      R(3,1) = 0.7745966692
      R(3,2) = 0.0
      R(3,3) = -R(3,1)
      CO(3,1) = 0.5555555556
      CO(3,2) = 0.8888888889
      CO(3,3) = CO(3,1)
      R(4,1) = 0.8611363116
      R(4,2) = 0.3399810436
      R(4,3) = -R(4,2)
      R(4,4) = -R(4,1)
      CO(4,1) = 0.3478548451
      CO(4,2) = 0.6521451549
      CO(4,3) = CO(4,2)
      CO(4,4) = CO(4,1)
      R(5,1) = 0.9061798459
      R(5,2) = 0.5384693101
      R(5,3) = 0.0
      R(5,4) = -R(5,2)
      R(5,5) = -R(5,1)
      CO(5,1) = 0.2369268850
      CO(5,2) = 0.4786286705
      CO(5,3) = 0.5688888889
      CO(5,4) = CO(5,2)
      CO(5,5) = CO(5,1)
C        STEP 1
      H1 = ( B - A ) / 2.0
      H2 = ( B + A ) / 2.0
C     USE AJ INSTEAD OF J
      AJ = 0.0
C       STEP 2
         DO 20 I = 1,M
C              STEP 3
               X = H1 * R(M,I) + H2
               JX = 0.0
               C1 = C( X )
               D1 = D( X )
               K1 = ( D1 - C1 ) / 2.0
               K2 = ( D1 + C1 ) / 2.0
C              STEP 4
               DO 30 J = 1,N
                     Y = K1 * R(N,J) + K2
                     Q = F( X, Y )
                     JX = JX + CO(N,J) * Q
30             CONTINUE
C              STEP 5
               AJ = AJ + CO(M,I) * K1 * JX
20       CONTINUE
C        STEP 6
         AJ = AJ * H1
C        STEP 7
C     OUTPUT
      WRITE(6,1) A,B
      WRITE(6,2) AJ
      WRITE(6,3) N,M
400   CLOSE(UNIT=5)
      CLOSE(UNIT=6)
      STOP
1     FORMAT(1X,'The double integral of F from ',E15.8,
     * ' to ',E15.8,' is')
2     FORMAT(1X,E15.8)
3     FORMAT(1X,'obtained with N = ',I3,' and M = ',I3)
      END

une erreur à la compilation apparait:

Code :

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

1
2
run-time error M6201:MATH
exp :domain error

svp aidez moi pour connaitre la cause de cette erreur.
merci.

**Sylvain Bergeron** · 09/07/2009, 17h22

Ce n'est pas une erreur de compilation, mais d'exécution (run-time error).

Ceci dit, le message indique « exp :domain error ». Ça veut dire que tu passes un mauvais argument à la fonction exp(). Erreur de domaine veut dire que la valeur de l'argument est en dehors de la plage pour laquelle la fonction est définie (ex: log(-2.)). Pour la fonction exp(), je dirais que tu l'appelles avec un argument trop grand qui produirait une valeur non représentable, ou avec un NaN (not a number), résultant d'une division par zéro ?

Comme il n'y a qu'un appel à exp() dans ton code (à la ligne 136, via la fonction interne F()), ça devrait être assez facile de trouver le problème en ajoutant un write(*,*) x,y juste avant la ligne 136...

**suzanne1307** · 10/07/2009, 00h48

Bonsoir,
le problème est que dans cette fonction on fait la double intégrale de la fonction exp(xz/yz) par rapport à xz et yz, xz allant de a à b et yz allant de x**2 à x**3 donc on n'est pas supposé connaitre les valeurs de xz et yz puisqu'elles sont des variables donc la write(*,*)xz,yz ne donnera aucun résultat.
est ce qu'il existe une autre solution pour connaitre l'origine de cette erreur?
merci

**FR119492** · 10/07/2009, 14h41

Salut!
Ton problème est formulé de manière assez incompréhensible. Pour qu'on y voie plus clair, je l'ai reformulé:

Dans le plan x,y, on a un domaine limité par des segments des droites x=a et x=b, ainsi que par des segments des courbes y=x**2 et y=x**3. On cherche l'intégrale de exp(x/y) sur ce domaine.

Je pense que, plutôt que d'utiliser une fonction toute faite, tu devrais la programmer toi-même.
Jean-Marc Blanc

**Sylvain Bergeron** · 10/07/2009, 19h33

Envoyé par suzanne1307

Bonsoir,
le problème est que dans cette fonction on fait la double intégrale de la fonction exp(xz/yz) par rapport à xz et yz, xz allant de a à b et yz allant de x**2 à x**3 donc on n'est pas supposé connaitre les valeurs de xz et yz puisqu'elles sont des variables donc la write(*,*)xz,yz ne donnera aucun résultat.
est ce qu'il existe une autre solution pour connaitre l'origine de cette erreur?
merci

Tu n'as pas bien lu ce que j'ai écrit :

Le seul appel que tu fais à la fonction exp() se trouve à la ligne 43
Ce qui se trouve à la ligne 43 est un « statement function ». Ça défini la fonction F(XZ,YZ) comme étant EXP(YZ/XZ)
À la ligne 136, ton programme appelle f(x,y)
Si tu appelles exp() va un mauvais argument, tu connais comment est établi cet argument à la ligne 136
Tu dois donc ajouter un write (*,*) juste avant la ligne 136 pour imprimer les arguments de F(), doit x et y (et non xz et yz)

double intégrale fortran

Fortran

Discussions similaires

Partager

Partager