| Possible crippling bug in NR_(SIMP,ROMB) [message #4121] |
Tue, 02 May 1995 00:00 |
mhaffner
Messages: 2
Registered: February 1995
|
Junior Member |
|
|
I get erroneous results when using the NR_SIMP and NR_ROMB routines
recursively that I belive traces to the routines themselves. Here is one
of the simplest routines I could think of to illustrate the problem (the
idea for this recursive calling comes from NR, section 4.6):
=========
FUNCTION myfunc, x
return, 2.0*x
END
FUNCTION innerint, y
return, nr_qromb('myfunc', 0.0, 5.0)
END
PRO simpleint
print, nr_qromb('innerint', 0.0, 2.0)
END
=========
Now, executing a statement like 'print, nr_qromb('myfunc',0.0,5.0)'
correctly returns 25.000. However when adding the extra integral from 0-2
over y as shown above in 'simpleint', IDL returns 37.500 when the correct
answer should be 50. I have run this exact calculation through the C
versions of the NR routines and they correctly report 50.000. I have
checked this on SGI's, Alpha's, and DecStations with the same results.
It seems the recursive calling is mangling the results. (Strangely, using
the y range of 0-3 gives 50!) If I replace the return statment in
'innerint' with 'return, 25.0' I (of course) get 50. As listed above,
'innerint' should be returning 25.000 every call.
Of course the integral above doesn't need to be calculated this way, but I
have some that I need to do using the recursive trick. Has anyone
encountered this before and found a workaround? Am I missing something
obvious about recursive function calling in IDL?
Thanks.
mh
--
Matt Haffner mhaffner@fosters.astro.wisc.edu
Dept. of Astronomy University of Wisconsin, Madison
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
