integration.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: integration.f90 #
5!# #
6!# Copyright (C) 2006-2014 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
32!----------------------------------------------------------------------------!
33MODULE integration
34 IMPLICIT NONE
35 !--------------------------------------------------------------------------!
36 PRIVATE
37#ifndef GAUSSN
38#define GAUSSN (15)
39
41#endif
42#define GAUSSN2 INT(GAUSSN / 2.0)
43
45 INTEGER, PARAMETER :: maxrec = 100
46 ! exclude definition of abscissas and weights from documentation
47 ! because it confuses doxygen
49#if GAUSSN == 2
50 REAL, DIMENSION(GAUSSN), PARAMETER :: &
51 gaussxk = (/ -0.577350269189626, 0.577350269189626 /)
52 REAL, DIMENSION(GAUSSN), PARAMETER :: &
53 gaussak = (/ 1.000000000000000, 1.000000000000000 /)
54#elif GAUSSN == 4
55 REAL, DIMENSION(GAUSSN), PARAMETER :: &
56 gaussxk = (/ -0.861136311594053, -0.339981043584856, &
57 0.339981043584856, 0.861136311594053 /)
58 REAL, DIMENSION(GAUSSN), PARAMETER :: &
59 gaussak = (/ 0.347854845137454, 0.652145154862546, &
60 0.652145154862546, 0.347854845137454 /)
61#elif GAUSSN == 8
62 REAL, DIMENSION(GAUSSN), PARAMETER :: &
63 gaussxk = (/ -0.960289856497536, -0.796666477413627, &
64 -0.525532409916329, -0.183434642495650, &
65 0.183434642495650, 0.525532409916329, &
66 0.796666477413627, 0.960289856497536 /)
67 REAL, DIMENSION(GAUSSN), PARAMETER :: &
68 gaussak = (/ 0.101228536290376, 0.222381034453374, &
69 0.313706645877887, 0.362683783378362, &
70 0.362683783378362, 0.313706645877887, &
71 0.222381034453374, 0.101228536290376 /)
72#elif GAUSSN == 12
73 REAL, DIMENSION(GAUSSN), PARAMETER :: &
74 gaussxk = (/ -0.9815606342467193, -0.9041172563704749, &
75 -0.7699026741943047, -0.5873179542866174, &
76 -0.3678314989981802, -0.1252334085114689, &
77 0.1252334085114689, 0.3678314989981802, &
78 0.5873179542866174, 0.7699026741943047, &
79 0.9041172563704749, 0.9815606342467193 /)
80 REAL, DIMENSION(GAUSSN), PARAMETER :: &
81 gaussak = (/ 0.047175336386512, 0.106939325995318, &
82 0.160078328543346, 0.203167426723066, &
83 0.233492536538355, 0.249147045813403, &
84 0.249147045813403, 0.233492536538355, &
85 0.203167426723066, 0.160078328543346, &
86 0.106939325995318, 0.04717533638651 /)
87#elif GAUSSN == 15
88 REAL, DIMENSION(GAUSSN), PARAMETER :: &
89 gaussxk = (/ -0.9879925180204854, -0.9372733924007059, &
90 -0.8482065834104272, -0.7244177313601700, &
91 -0.5709721726085388, -0.3941513470775634, &
92 -0.2011940939974345, &
93 0.0, &
94 0.2011940939974345, &
95 0.3941513470775634, 0.5709721726085388, &
96 0.7244177313601700, 0.8482065834104272, &
97 0.9372733924007059, 0.9879925180204854 /)
98 REAL, DIMENSION(GAUSSN), PARAMETER :: &
99 gaussak = (/ 0.030753241996117, 0.070366047488108, &
100 0.107159220467172, 0.139570677926154, &
101 0.166269205816994, 0.186161000015562, &
102 0.198431485327112, &
103 0.202578241925561, &
104 0.198431485327112, &
105 0.186161000015562, 0.166269205816994, &
106 0.139570677926154, 0.107159220467172, &
107 0.070366047488108, 0.03075324199612 /)
108#else
109# error Wrong GAUSSN number in numtools/integration.f90
110#endif
111
112 REAL, DIMENSION(3,MAXREC) :: stack
113 !--------------------------------------------------------------------------!
114 PUBLIC integrate
115 !--------------------------------------------------------------------------!
116
117CONTAINS
118
124 FUNCTION integrate(fkt,xl,xr,eps,plist,method) RESULT(integral)
125 IMPLICIT NONE
126 !----------------------------------------------------------------------!
127 INTERFACE
128 FUNCTION fkt(x,plist) RESULT(fx)
129 IMPLICIT NONE
130 REAL, INTENT(IN) :: x
131 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
132 REAL :: fx
133 END FUNCTION fkt
134 END INTERFACE
136 !----------------------------------------------------------------------!
137 REAL :: xl
138 REAL :: xr
139 REAL :: eps
140 REAL, DIMENSION(:), OPTIONAL :: &
141 plist
142 INTEGER, OPTIONAL :: method
143 REAL :: integral
144 !----------------------------------------------------------------------!
145 INTEGER :: meth
146 REAL :: tmps,err
147 !----------------------------------------------------------------------!
148 INTENT(IN) :: xl,xr,eps,method
149 INTENT(INOUT) :: plist
150 !----------------------------------------------------------------------!
151
152 IF (xl.EQ.xr) THEN
153 integral = 0.
154 RETURN
155 END IF
156
161 IF (PRESENT(method).AND.((method.GT.0).AND.method.LE.3)) THEN
162 meth = method
163 ELSE
164 ! default is recursive Gauss
165 meth = 1
166 END IF
167
168 ! tolerance for integration
169 err = abs(xr-xl) * eps
170
171 SELECT CASE(meth)
172 CASE(1)
173 ! recursive Gauss integration
174 tmps = 0.
175 IF (PRESENT(plist)) THEN
176 tmps = qgauss1d(fkt,min(xl,xr),max(xl,xr),plist)
177 integral = sign(1.0,xr-xl)*qadaptive1d_recursive(fkt,tmps,&
178 min(xl,xr),max(xl,xr),err,plist)
179 ELSE
180 tmps = qgauss1d(fkt,min(xl,xr),max(xl,xr))
181 integral = sign(1.0,xr-xl)*qadaptive1d_recursive(fkt,tmps,&
182 min(xl,xr),max(xl,xr),err)
183 END IF
184
185 CASE(2)
186 ! iterative Gauss integration
187 tmps = 0.
188 IF (PRESENT(plist)) THEN
189 tmps = qgauss1d(fkt,min(xl,xr),max(xl,xr),plist)
190 integral = sign(1.0,xr-xl)*qadaptive1d_iterative(fkt,tmps,&
191 min(xl,xr),max(xl,xr),err,plist)
192 ELSE
193 tmps = qgauss1d(fkt,min(xl,xr),max(xl,xr))
194 integral = sign(1.0,xr-xl)*qadaptive1d_iterative(fkt,tmps,&
195 min(xl,xr),max(xl,xr),err)
196 END IF
197
198 CASE(3)
199 ! Romberg integration
200 IF (PRESENT(plist)) THEN
201 integral = sign(1.0,xr-xl)*qromberg(fkt,min(xl,xr),max(xl,xr),err,plist)
202 ELSE
203 integral = sign(1.0,xr-xl)*qromberg(fkt,min(xl,xr),max(xl,xr),err)
204 END IF
205
206 CASE DEFAULT
207 print *, "ERROR in integrate: select one of 1,2 for method"
208 stop
209 END SELECT
210
211 END FUNCTION integrate
212
213
214 !------------------------------------------------------------------------!
215 !
216 ! 1D Integration with adaptive Gauss
217 !
218 !------------------------------------------------------------------------!
219 RECURSIVE FUNCTION qadaptive1d_recursive(fkt, oldS, xl, xr, tol, plist) RESULT(S)
220 IMPLICIT NONE
221 !----------------------------------------------------------------------!
222 INTERFACE
223 FUNCTION fkt(x,plist) RESULT(fx)
224 IMPLICIT NONE
225 REAL, INTENT(IN) :: x
226 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
227 REAL :: fx
228 END FUNCTION fkt
229 END INTERFACE
230 !----------------------------------------------------------------------!
231 REAL :: olds, xl, xr, tol
232 REAL, DIMENSION(:), OPTIONAL :: plist
233 REAL :: s
234 !----------------------------------------------------------------------!
235 REAL :: m, res, resl, resr, qerr
236 !----------------------------------------------------------------------!
237 INTENT(IN) :: olds, xl, xr, tol
238 INTENT(INOUT) :: plist
239 !----------------------------------------------------------------------!
240
241 m = 0.5 * (xr+xl)
242 resl = qgauss1d(fkt, xl, m, plist) ! left solution
243 resr = qgauss1d(fkt, m, xr, plist) ! right solution
244 res = resl + resr
245
246 ! error estimate
247 qerr = (olds - res) / (4**gaussn - 1.0)
248
249 IF (abs(qerr).GE.tol) THEN
250 ! recursion
251 s = qadaptive1d_recursive(fkt, resl, xl, m, 0.5*tol, plist) &
252 + qadaptive1d_recursive(fkt, resr, m, xr, 0.5*tol, plist)
253 ELSE
254 s = res + qerr
255 END IF
256 END FUNCTION qadaptive1d_recursive
257
258
259 FUNCTION qadaptive1d_iterative(fkt, oldS, xl, xr, tol, plist) RESULT(S)
260 IMPLICIT NONE
261 !----------------------------------------------------------------------!
262 INTERFACE
263 FUNCTION fkt(x,plist) RESULT(fx)
264 IMPLICIT NONE
265 REAL, INTENT(IN) :: x
266 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
267 REAL :: fx
268 END FUNCTION fkt
269 END INTERFACE
270 !----------------------------------------------------------------------!
271 REAL :: olds, xl, xr, tol
272 REAL, DIMENSION(:), OPTIONAL :: plist
273 REAL :: s
274 !----------------------------------------------------------------------!
275 INTEGER :: i,sptr
276 REAL :: resl, resr, reso, res, qerr
277 REAL :: l,r,m,err
278 !----------------------------------------------------------------------!
279 INTENT(IN) :: olds, tol, xl, xr
280 INTENT(INOUT) :: plist
281 !----------------------------------------------------------------------!
282
283 ! copy basic variables
284 l = xl
285 r = xr
286 reso = olds
287 err = tol
288
289 ! initialize result and stack pointer
290 s = 0.
291 sptr = 1
292
293 ! main loop
294 DO i=1,maxrec
295 m = 0.5 * (r+l) ! devide interval
296 resl = qgauss1d(fkt, l, m, plist) ! integrate over [xl..m]
297 resr = qgauss1d(fkt, m, r, plist) ! integrate over [m..xr]
298 res = resl + resr ! add two results
299 qerr = (reso - res) / (4**gaussn - 1.0) ! error estimate
300 ! check error
301 IF (abs(qerr) < err) THEN
302 ! store result
303 s = s + res + qerr
304 IF (sptr.GT.1) THEN
305 ! pop values from stack
306 sptr = sptr-1
307 l = stack(1,sptr)
308 r = stack(2,sptr)
309 reso = stack(3,sptr)
310 ELSE
311 ! finished
312 EXIT
313 END IF
314 ELSE
315 ! push right values on stack
316 stack(1,sptr) = m
317 stack(2,sptr) = r
318 stack(3,sptr) = resr
319 sptr = sptr+1
320 ! continue with left
321 r = m
322 err = 0.5*err
323 END IF
324 END DO
325
326 IF (i.GE.maxrec) THEN
327 print *, "ERROR in qadaptive1D_iterative: max recursion reached"
328 stop
329 END IF
330
331 END FUNCTION qadaptive1d_iterative
332
333
334 !------------------------------------------------------------------------!
335 !
336 ! basic 1D Gauss quadrature function
337 !
338 !------------------------------------------------------------------------!
339 FUNCTION qgauss1d(fkt,xl,xr,plist) RESULT(Sgauss)
340 IMPLICIT NONE
341 !----------------------------------------------------------------------!
342 INTERFACE
343 FUNCTION fkt(x,plist) RESULT(fx)
344 IMPLICIT NONE
345 REAL, INTENT(IN) :: x
346 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
347 REAL :: fx
348 END FUNCTION fkt
349 END INTERFACE
350 !----------------------------------------------------------------------!
351 REAL, INTENT(IN) :: xl, xr
352 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist ! additional parameters
353 REAL :: sgauss
354 REAL :: m, h, t
355 INTEGER :: j
356 !----------------------------------------------------------------------!
357
358 m = 0.5 * (xr+xl)
359 h = 0.5 * (xr-xl)
360 sgauss = 0.0
361
362 DO j=1,gaussn2
363 t = h * gaussxk(j)
364 sgauss = sgauss + gaussak(j) * (fkt(m-t,plist)+fkt(m+t,plist))
365 END DO
366
367 IF (mod(gaussn,2) /= 0) THEN
368 sgauss = sgauss + gaussak(gaussn2+1) * fkt(m,plist)
369 END IF
370 sgauss = h * sgauss
371 END FUNCTION qgauss1d
372
373
374 !------------------------------------------------------------------------!
375 !
376 ! 1D Integration with Romberg algorithm
377 !
378 !------------------------------------------------------------------------!
379 FUNCTION qromberg(fkt, xl, xr, tol, plist) RESULT (Srom)
380 IMPLICIT NONE
381 !----------------------------------------------------------------------!
382 INTERFACE
383 FUNCTION fkt(x, plist) RESULT(fx)
384 IMPLICIT NONE
385 REAL, INTENT(IN) :: x
386 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
387 REAL :: fx
388 END FUNCTION fkt
389 END INTERFACE
390 !----------------------------------------------------------------------!
391 REAL, INTENT(IN) :: xl, xr
392 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist ! additional parameters
393 REAL :: srom
394 INTEGER, PARAMETER :: jmax = 10 ! max refinement
395 REAL, DIMENSION(jmax,jmax) :: res
396 REAL :: h
397 REAL :: tol
398 INTEGER :: i, j, n
399 !----------------------------------------------------------------------!
400
401 n = 2
402 ! step size
403 h = (xr - xl) / n
404 res(:,:) = 0.0
405 ! sum of trapezoidal rule at the boundaries
406 res(1,1) = 0.5 * (fkt(xl,plist) + fkt(xr,plist))
407 ! and in between
408 DO i=1,n-1
409 res(1,1) = res(1,1) + fkt(xl+i*h,plist)
410 END DO
411 res(1,1) = res(1,1) * h
412
413 DO j=2,jmax
414 ! half step size
415 h = 0.5 * h
416 ! double grid points
417 n = 2 * n
418 ! sum up contributions due to the new points
419 DO i=1,n-1,2
420 res(j,1) = res(j,1) + fkt(xl+i*h,plist)
421 END DO
422 res(j,1) = 0.5 * res(j-1,1) + h*res(j,1)
423 ! extrapolation
424 DO i=1,j-1
425 res(j,i+1) = res(j,i) + (res(j,i) - res(j-1,i)) / (4**i - 1)
426 END DO
427 ! check convergence
428 IF (abs(res(j,j)-res(j,j-1))<tol) EXIT
429 END DO
430
431 srom = res(j-1,j-1)
432 END FUNCTION qromberg
433
434END MODULE integration
integration
Numerical integration.
Definition: integration.f90:33
integration::qgauss1d
real function qgauss1d(fkt, xl, xr, plist)
Definition: integration.f90:340
integration::maxrec
integer, parameter maxrec
maximal depth for recursion
Definition: integration.f90:45
integration::integrate
real function, public integrate(fkt, xl, xr, eps, plist, method)
Numerical integration function.
Definition: integration.f90:125
integration::qromberg
real function qromberg(fkt, xl, xr, tol, plist)
Definition: integration.f90:380
integration::stack
real, dimension(3, maxrec) stack
stack for iterative algorithm
Definition: integration.f90:112
integration::qadaptive1d_iterative
real function qadaptive1d_iterative(fkt, oldS, xl, xr, tol, plist)
Definition: integration.f90:260
integration::qadaptive1d_recursive
recursive real function qadaptive1d_recursive(fkt, oldS, xl, xr, tol, plist)
Definition: integration.f90:220
Generated by   doxygen 1.9.4