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
dqk21.f
Go to the documentation of this file.
1  SUBROUTINE dqk21(F,A,B,RESULT,ABSERR,RESABS,RESASC,IERR)
2 C***BEGIN PROLOGUE DQK21
3 C***DATE WRITTEN 800101 (YYMMDD)
4 C***REVISION DATE 830518 (YYMMDD)
5 C***CATEGORY NO. H2A1A2
6 C***KEYWORDS 21-POINT GAUSS-KRONROD RULES
7 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
8 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
9 C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
10 C ESTIMATE
11 C J = INTEGRAL OF ABS(F) OVER (A,B)
12 C***DESCRIPTION
13 C
14 C INTEGRATION RULES
15 C STANDARD FORTRAN SUBROUTINE
16 C DOUBLE PRECISION VERSION
17 C
18 C PARAMETERS
19 C ON ENTRY
20 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
21 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
22 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
23 C
24 C A - DOUBLE PRECISION
25 C LOWER LIMIT OF INTEGRATION
26 C
27 C B - DOUBLE PRECISION
28 C UPPER LIMIT OF INTEGRATION
29 C
30 C ON RETURN
31 C RESULT - DOUBLE PRECISION
32 C APPROXIMATION TO THE INTEGRAL I
33 C RESULT IS COMPUTED BY APPLYING THE 21-POINT
34 C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION
35 C OF ABSCISSAE TO THE 10-POINT GAUSS RULE (RESG).
36 C
37 C ABSERR - DOUBLE PRECISION
38 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
39 C WHICH SHOULD NOT EXCEED ABS(I-RESULT)
40 C
41 C RESABS - DOUBLE PRECISION
42 C APPROXIMATION TO THE INTEGRAL J
43 C
44 C RESASC - DOUBLE PRECISION
45 C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
46 C OVER (A,B)
47 C
48 C***REFERENCES (NONE)
49 C***ROUTINES CALLED D1MACH
50 C***END PROLOGUE DQK21
51 C
52  DOUBLE PRECISION a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
53  * d1mach,epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
54  * resg,resk,reskh,result,uflow,wg,wgk,xgk
55  INTEGER j,jtw,jtwm1
56  EXTERNAL f
57 C
58  dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
59 C
60 C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
61 C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
62 C CORRESPONDING WEIGHTS ARE GIVEN.
63 C
64 C XGK - ABSCISSAE OF THE 21-POINT KRONROD RULE
65 C XGK(2), XGK(4), ... ABSCISSAE OF THE 10-POINT
66 C GAUSS RULE
67 C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
68 C ADDED TO THE 10-POINT GAUSS RULE
69 C
70 C WGK - WEIGHTS OF THE 21-POINT KRONROD RULE
71 C
72 C WG - WEIGHTS OF THE 10-POINT GAUSS RULE
73 C
74 C
75 C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS
76 C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
77 C BELL LABS, NOV. 1981.
78 C
79  DATA wg( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
80  DATA wg( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
81  DATA wg( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
82  DATA wg( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
83  DATA wg( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
84 C
85  DATA xgk( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
86  DATA xgk( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
87  DATA xgk( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
88  DATA xgk( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
89  DATA xgk( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
90  DATA xgk( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
91  DATA xgk( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
92  DATA xgk( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
93  DATA xgk( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
94  DATA xgk( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
95  DATA xgk( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
96 C
97  DATA wgk( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
98  DATA wgk( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
99  DATA wgk( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
100  DATA wgk( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
101  DATA wgk( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
102  DATA wgk( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
103  DATA wgk( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
104  DATA wgk( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
105  DATA wgk( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
106  DATA wgk( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
107  DATA wgk( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
108 C
109 C
110 C LIST OF MAJOR VARIABLES
111 C -----------------------
112 C
113 C CENTR - MID POINT OF THE INTERVAL
114 C HLGTH - HALF-LENGTH OF THE INTERVAL
115 C ABSC - ABSCISSA
116 C FVAL* - FUNCTION VALUE
117 C RESG - RESULT OF THE 10-POINT GAUSS FORMULA
118 C RESK - RESULT OF THE 21-POINT KRONROD FORMULA
119 C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
120 C I.E. TO I/(B-A)
121 C
122 C
123 C MACHINE DEPENDENT CONSTANTS
124 C ---------------------------
125 C
126 C EPMACH IS THE LARGEST RELATIVE SPACING.
127 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
128 C
129 C***FIRST EXECUTABLE STATEMENT DQK21
130  epmach = d1mach(4)
131  uflow = d1mach(1)
132 C
133  centr = 0.5d+00*(a+b)
134  hlgth = 0.5d+00*(b-a)
135  dhlgth = dabs(hlgth)
136 C
137 C COMPUTE THE 21-POINT KRONROD APPROXIMATION TO
138 C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
139 C
140  resg = 0.0d+00
141  ierr = 0
142  CALL f(centr,ierr,fc)
143  IF (ierr .LT. 0) RETURN
144  resk = wgk(11)*fc
145  resabs = dabs(resk)
146  DO 10 j=1,5
147  jtw = 2*j
148  absc = hlgth*xgk(jtw)
149  CALL f(centr-absc,ierr,fval1)
150  IF (ierr .LT. 0) RETURN
151  CALL f(centr+absc,ierr,fval2)
152  IF (ierr .LT. 0) RETURN
153  fv1(jtw) = fval1
154  fv2(jtw) = fval2
155  fsum = fval1+fval2
156  resg = resg+wg(j)*fsum
157  resk = resk+wgk(jtw)*fsum
158  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
159  10 CONTINUE
160  DO 15 j = 1,5
161  jtwm1 = 2*j-1
162  absc = hlgth*xgk(jtwm1)
163  CALL f(centr-absc,ierr,fval1)
164  IF (ierr .LT. 0) RETURN
165  CALL f(centr+absc,ierr,fval2)
166  IF (ierr .LT. 0) RETURN
167  fv1(jtwm1) = fval1
168  fv2(jtwm1) = fval2
169  fsum = fval1+fval2
170  resk = resk+wgk(jtwm1)*fsum
171  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
172  15 CONTINUE
173  reskh = resk*0.5d+00
174  resasc = wgk(11)*dabs(fc-reskh)
175  DO 20 j=1,10
176  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
177  20 CONTINUE
178  result = resk*hlgth
179  resabs = resabs*dhlgth
180  resasc = resasc*dhlgth
181  abserr = dabs((resk-resg)*hlgth)
182  IF(resasc.NE.0.0d+00.AND.abserr.NE.0.0d+00)
183  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
184  IF(resabs.GT.uflow/(0.5d+02*epmach)) abserr = dmax1
185  * ((epmach*0.5d+02)*resabs,abserr)
186  RETURN
187  END
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2