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 |
Partager