| Bug (still) in NR_QROMB, QROMB [message #4981] |
Wed, 30 August 1995 00:00 |
mallozzi
Messages: 60
Registered: August 1994
|
Member |
|
|
; There appears to be a bug in the routines NR_QROMB and QROMB if you call
; the functions recursively. I included a test program below.
; I notified RSI several months ago about this, but apparently the error
; never got corrected. Be careful how you construct the calling sequences
; for these routines!
;
; If you are running IDL 4.0 and want to test the function QROMB (not
; available before 4.0) then uncomment the appropriate cases below.
;
; -Bob Mallozzi
;
;
;- - - - - - - - - - - - - - - - - - - CUT - - - - - - - - - - - - - - - - - -
; Here is a sample program I used to test NR_QROMB.
; It computes an integral whose integrand contains a second integral:
; I = int_0^2{ x^2 * [ int_0^1{y^3 * dy} ] * dx}
; = 0.66666
;
; NOTE THE RECURSIVE CALLING SEQUENCE, WHICH PRODUCES THE ERROR.
;
; Calls to SIMPSON produce the correct result.
; Calls to NR_QROMB, QROMB produce the correct result only if
; one of the calls (it doesn't matter which one) is DOUBLE precision.
; If both calls are DOUBLE or both calls are FLOAT, the result is
; not correct.
;
; Tested on following machines:
; !VERSION = { vax vms 3.5.1}
; !VERSION = { alpha vms 3.6}
; !VERSION = { alpha vms vms 4.0.1}
;
; I tested several combinations of the integral limits; this
; did not affect the results.
;
; Robert Mallozzi
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function func1, x
common data, inner, select
; Get the inner integral
ll = 0.0
ul = 1.0
case select of
1: inner = SIMPSON('func2', ll, ul, TOL = 0.001)
2: inner = SIMPSON('func2', DOUBLE(ll), DOUBLE(ul), TOL = 0.001)
3: inner = SIMPSON('func2', DOUBLE(ll), DOUBLE(ul), TOL = 0.001)
4: inner = NR_QROMB('func2', ll, ul, EPS = 0.001)
5: inner = NR_QROMB('func2', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
6: inner = NR_QROMB('func2', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
; 7: inner = QROMB('func2', ll, ul, EPS = 0.001)
; 8: inner = QROMB('func2', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
; 9: inner = QROMB('func2', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
endcase
return, (x^2) * inner
end
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function func2, y
return, y^3
end
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
common data, inner, select
choices = ['Select the integration routines: ', $
'SIMPSON single, SIMPSON single', $
'SIMPSON double, SIMPSON double', $
'SIMPSON single, SIMPSON double', $
'NR_QROMB single, NR_QROMB single', $
'NR_QROMB double, NR_QROMB double', $
'NR_QROMB single, NR_QROMB double']
; 'QROMB single, QROMB single', $
; 'QROMB double, QROMB double', $
; 'QROMB single, QROMB double']
select = WMENU(choices, title = 0, init = 1)
ll = 0.0
ul = 2.0
case select of
1: BEGIN
result = SIMPSON('func1', ll, ul, TOL = 0.001)
print, choices(1)
END
2: BEGIN
result = SIMPSON('func1', DOUBLE(ll), DOUBLE(ul), TOL = 0.001)
print, choices(3)
END
3: BEGIN
result = SIMPSON('func1', ll, ul, TOL = 0.001)
print, choices(2)
END
4: BEGIN
result = NR_QROMB('func1', ll, ul, EPS = 0.001)
print, choices(4)
END
5: BEGIN
result = NR_QROMB('func1', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
print, choices(5)
END
6: BEGIN
result = NR_QROMB('func1', ll, ul, EPS = 0.001)
print, choices(6)
END
; 7: BEGIN
; result = QROMB('func1', ll, ul, EPS = 0.001)
; print, choices(7)
; END
; 8: BEGIN
; result = QROMB('func1', DOUBLE(ll), DOUBLE(ul), EPS = 0.001)
; print, choices(8)
; END
; 9: BEGIN
; result = QROMB('func1', ll, ul, EPS = 0.001)
; print, choices(9)
; END
endcase
print, ' Numerical result is ', result
print, ' Analytical result is ', (1.0 / 3) * (ul^3 - ll^3) * inner
print
END
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
