GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dqk15i.f
Go to the documentation of this file.
1  SUBROUTINE dqk15i(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC,
2  1 ierr)
3 C***BEGIN PROLOGUE DQK15I
4 C***DATE WRITTEN 800101 (YYMMDD)
5 C***REVISION DATE 830518 (YYMMDD)
6 C***CATEGORY NO. H2A3A2,H2A4A2
7 C***KEYWORDS 15-POINT TRANSFORMED GAUSS-KRONROD RULES
8 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
9 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
10 C***PURPOSE THE ORIGINAL (INFINITE INTEGRATION RANGE IS MAPPED
11 C ONTO THE INTERVAL (0,1) AND (A,B) IS A PART OF (0,1).
12 C IT IS THE PURPOSE TO COMPUTE
13 C I = INTEGRAL OF TRANSFORMED INTEGRAND OVER (A,B),
14 C J = INTEGRAL OF ABS(TRANSFORMED INTEGRAND) OVER (A,B).
15 C***DESCRIPTION
16 C
17 C INTEGRATION RULE
18 C STANDARD FORTRAN SUBROUTINE
19 C DOUBLE PRECISION VERSION
20 C
21 C PARAMETERS
22 C ON ENTRY
23 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
24 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
25 C DECLARED E X T E R N A L IN THE CALLING PROGRAM.
26 C
27 C BOUN - DOUBLE PRECISION
28 C FINITE BOUND OF ORIGINAL INTEGRATION
29 C RANGE (SET TO ZERO IF INF = +2)
30 C
31 C INF - INTEGER
32 C IF INF = -1, THE ORIGINAL INTERVAL IS
33 C (-INFINITY,BOUND),
34 C IF INF = +1, THE ORIGINAL INTERVAL IS
35 C (BOUND,+INFINITY),
36 C IF INF = +2, THE ORIGINAL INTERVAL IS
37 C (-INFINITY,+INFINITY) AND
38 C THE INTEGRAL IS COMPUTED AS THE SUM OF TWO
39 C INTEGRALS, ONE OVER (-INFINITY,0) AND ONE OVER
40 C (0,+INFINITY).
41 C
42 C A - DOUBLE PRECISION
43 C LOWER LIMIT FOR INTEGRATION OVER SUBRANGE
44 C OF (0,1)
45 C
46 C B - DOUBLE PRECISION
47 C UPPER LIMIT FOR INTEGRATION OVER SUBRANGE
48 C OF (0,1)
49 C
50 C ON RETURN
51 C RESULT - DOUBLE PRECISION
52 C APPROXIMATION TO THE INTEGRAL I
53 C RESULT IS COMPUTED BY APPLYING THE 15-POINT
54 C KRONROD RULE(RESK) OBTAINED BY OPTIMAL ADDITION
55 C OF ABSCISSAE TO THE 7-POINT GAUSS RULE(RESG).
56 C
57 C ABSERR - DOUBLE PRECISION
58 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
59 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
60 C
61 C RESABS - DOUBLE PRECISION
62 C APPROXIMATION TO THE INTEGRAL J
63 C
64 C RESASC - DOUBLE PRECISION
65 C APPROXIMATION TO THE INTEGRAL OF
66 C ABS((TRANSFORMED INTEGRAND)-I/(B-A)) OVER (A,B)
67 C
68 C***REFERENCES (NONE)
69 C***ROUTINES CALLED D1MACH
70 C***END PROLOGUE DQK15I
71 C
72  DOUBLE PRECISION a,absc,absc1,absc2,abserr,b,boun,centr,dabs,dinf,
73  * dmax1,dmin1,d1mach,epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,
74  * resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,uflow,wg,wgk,
75  * xgk,fvalt
76  INTEGER inf,j
77  EXTERNAL f
78 C
79  dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
80 C
81 C THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
82 C (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
83 C THEIR CORRESPONDING WEIGHTS ARE GIVEN.
84 C
85 C XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE
86 C XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
87 C GAUSS RULE
88 C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
89 C ADDED TO THE 7-POINT GAUSS RULE
90 C
91 C WGK - WEIGHTS OF THE 15-POINT KRONROD RULE
92 C
93 C WG - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
94 C TO THE ABSCISSAE XGK(2), XGK(4), ...
95 C WG(1), WG(3), ... ARE SET TO ZERO.
96 C
97  DATA wg(1) / 0.0d0 /
98  DATA wg(2) / 0.1294849661 6886969327 0611432679 082d0 /
99  DATA wg(3) / 0.0d0 /
100  DATA wg(4) / 0.2797053914 8927666790 1467771423 780d0 /
101  DATA wg(5) / 0.0d0 /
102  DATA wg(6) / 0.3818300505 0511894495 0369775488 975d0 /
103  DATA wg(7) / 0.0d0 /
104  DATA wg(8) / 0.4179591836 7346938775 5102040816 327d0 /
105 C
106  DATA xgk(1) / 0.9914553711 2081263920 6854697526 329d0 /
107  DATA xgk(2) / 0.9491079123 4275852452 6189684047 851d0 /
108  DATA xgk(3) / 0.8648644233 5976907278 9712788640 926d0 /
109  DATA xgk(4) / 0.7415311855 9939443986 3864773280 788d0 /
110  DATA xgk(5) / 0.5860872354 6769113029 4144838258 730d0 /
111  DATA xgk(6) / 0.4058451513 7739716690 6606412076 961d0 /
112  DATA xgk(7) / 0.2077849550 0789846760 0689403773 245d0 /
113  DATA xgk(8) / 0.0000000000 0000000000 0000000000 000d0 /
114 C
115  DATA wgk(1) / 0.0229353220 1052922496 3732008058 970d0 /
116  DATA wgk(2) / 0.0630920926 2997855329 0700663189 204d0 /
117  DATA wgk(3) / 0.1047900103 2225018383 9876322541 518d0 /
118  DATA wgk(4) / 0.1406532597 1552591874 5189590510 238d0 /
119  DATA wgk(5) / 0.1690047266 3926790282 6583426598 550d0 /
120  DATA wgk(6) / 0.1903505780 6478540991 3256402421 014d0 /
121  DATA wgk(7) / 0.2044329400 7529889241 4161999234 649d0 /
122  DATA wgk(8) / 0.2094821410 8472782801 2999174891 714d0 /
123 C
124 C
125 C LIST OF MAJOR VARIABLES
126 C -----------------------
127 C
128 C CENTR - MID POINT OF THE INTERVAL
129 C HLGTH - HALF-LENGTH OF THE INTERVAL
130 C ABSC* - ABSCISSA
131 C TABSC* - TRANSFORMED ABSCISSA
132 C FVAL* - FUNCTION VALUE
133 C RESG - RESULT OF THE 7-POINT GAUSS FORMULA
134 C RESK - RESULT OF THE 15-POINT KRONROD FORMULA
135 C RESKH - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
136 C INTEGRAND OVER (A,B), I.E. TO I/(B-A)
137 C
138 C MACHINE DEPENDENT CONSTANTS
139 C ---------------------------
140 C
141 C EPMACH IS THE LARGEST RELATIVE SPACING.
142 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
143 C
144 C***FIRST EXECUTABLE STATEMENT DQK15I
145  epmach = d1mach(4)
146  uflow = d1mach(1)
147  dinf = min0(1,inf)
148 C
149  centr = 0.5d+00*(a+b)
150  hlgth = 0.5d+00*(b-a)
151  tabsc1 = boun+dinf*(0.1d+01-centr)/centr
152  ierr = 0
153  CALL f(tabsc1,ierr,fval1)
154  IF (ierr .LT. 0) RETURN
155  IF(inf.EQ.2) THEN
156  CALL f(-tabsc1,ierr,fvalt)
157  IF (ierr .LT. 0) RETURN
158  fval1 = fval1+fvalt
159  ENDIF
160  fc = (fval1/centr)/centr
161 C
162 C COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
163 C THE INTEGRAL, AND ESTIMATE THE ERROR.
164 C
165  resg = wg(8)*fc
166  resk = wgk(8)*fc
167  resabs = dabs(resk)
168  DO 10 j=1,7
169  absc = hlgth*xgk(j)
170  absc1 = centr-absc
171  absc2 = centr+absc
172  tabsc1 = boun+dinf*(0.1d+01-absc1)/absc1
173  tabsc2 = boun+dinf*(0.1d+01-absc2)/absc2
174  CALL f(tabsc1,ierr,fval1)
175  IF (ierr .LT. 0) RETURN
176  CALL f(tabsc2,ierr,fval2)
177  IF (ierr .LT. 0) RETURN
178  IF(inf.EQ.2) THEN
179  CALL f(-tabsc1,ierr,fvalt)
180  IF (ierr .LT. 0) RETURN
181  fval1 = fval1+fvalt
182  ENDIF
183  IF(inf.EQ.2) THEN
184  CALL f(-tabsc2,ierr,fvalt)
185  IF (ierr .LT. 0) RETURN
186  fval2 = fval2+fvalt
187  ENDIF
188  fval1 = (fval1/absc1)/absc1
189  fval2 = (fval2/absc2)/absc2
190  fv1(j) = fval1
191  fv2(j) = fval2
192  fsum = fval1+fval2
193  resg = resg+wg(j)*fsum
194  resk = resk+wgk(j)*fsum
195  resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2))
196  10 CONTINUE
197  reskh = resk*0.5d+00
198  resasc = wgk(8)*dabs(fc-reskh)
199  DO 20 j=1,7
200  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
201  20 CONTINUE
202  result = resk*hlgth
203  resasc = resasc*hlgth
204  resabs = resabs*hlgth
205  abserr = dabs((resk-resg)*hlgth)
206  IF(resasc.NE.0.0d+00.AND.abserr.NE.0.d0) abserr = resasc*
207  * dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
208  IF(resabs.GT.uflow/(0.5d+02*epmach)) abserr = dmax1
209  * ((epmach*0.5d+02)*resabs,abserr)
210  RETURN
211  END
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2