Fun with Int_tabulated.pro [message #68497] |
Fri, 30 October 2009 10:20  |
wlandsman
Messages: 743 Registered: June 2000
|
Senior Member |
|
|
I ran into a couple of gotchas with int_tabulated.pro this
morning, It performs numerical integration for the simple case:
IDL> x = findgen(5)
IDL> y = x^2
IDL> print,int_tabulated(x,y)
21.3333
But if one has flipped both vectors (say while converting wavelength
to frequency), then -- even though a plot of the two vectors looks
exactly the same -- int_tabulated gives a different answer.
IDL> x = reverse(x) & y = reverse(y)
IDL> print,int_tabulated(x,y)
42.3619
A closer look at the int_tabulated documentation tells us that if X if
not monotonic increasing, then one needs to set the /SORT keyword.
So now we follow the documentation:
IDL> print,int_tabulated(x,y,/sort)
21.3333
and get the right answer again. But as the documentation also "sort"
of warns
you, setting the /SORT keyword reorders the input X and Y variables.
This is
probably not a problem unless one has modified the Y vector on the fly
by
some factor f
f = findgen(5) + 0.5
print,int_tabulated(x, y*f,/sort)
Upon output, the X vector has been sorted, but the Y vector has not,
since it
was passed as part of a temporary variable y*f. So one has lost
the one to
one correspondence between X and Y, screwing up all subsequent
processing.
--Wayne
|
|
|
Re: Fun with Int_tabulated.pro [message #68534 is a reply to message #68497] |
Wed, 04 November 2009 09:50  |
Chris[6]
Messages: 84 Registered: July 2008
|
Member |
|
|
> Wow. That's one rule my Numerical Analysis prof drilled into our heads
> (too) many years ago: never, ever, ever, EVER modify the inputs to a
> function call. Create a temporary variable inside the function and do
> your modifications on that.
>
> Ed
I agree with this idea - any chance int_tabulated will be changed to
protect the input variables?
Chris
|
|
|
Re: Fun with Int_tabulated.pro [message #68536 is a reply to message #68497] |
Wed, 04 November 2009 08:47  |
edward.s.meinel@aero.
Messages: 52 Registered: February 2005
|
Member |
|
|
On Oct 30, 12:20 pm, wlandsman <wlands...@gmail.com> wrote:
> and get the right answer again. But as the documentation also "sort"
> of warns
> you, setting the /SORT keyword reorders the input X and Y variables.
>
> --Wayne
Wow. That's one rule my Numerical Analysis prof drilled into our heads
(too) many years ago: never, ever, ever, EVER modify the inputs to a
function call. Create a temporary variable inside the function and do
your modifications on that.
Ed
|
|
|
Re: Fun with Int_tabulated.pro [message #68546 is a reply to message #68497] |
Tue, 03 November 2009 16:18  |
Paul Van Delst[1]
Messages: 1157 Registered: April 2002
|
Senior Member |
|
|
Hello,
I've been struggling with int_tabulated recently but for a different reason: it always
performs spline interpolation. If you are integrating high-order Gaussian-like functions
with an insufficient point density (e.g. some satellite sensor measured spectral response
functions), it will always produce terrible results due to Gibbs-type phenomena. In
general that's not unexpected, but a bugger of a problem to track down in automated
day-to-day processing.
I'm writing my own integrator to allow for different interpolation schemes. It handles the
problem you're having too by simply copying the input to temporaries, and then doing all
the sorting/uniq'ing/etc.
Anyway.... apologies for hijacking the thread. I just wanted to mini-vent. :o)
cheers,
paulv
wlandsman wrote:
> I ran into a couple of gotchas with int_tabulated.pro this
> morning, It performs numerical integration for the simple case:
>
> IDL> x = findgen(5)
> IDL> y = x^2
> IDL> print,int_tabulated(x,y)
> 21.3333
>
> But if one has flipped both vectors (say while converting wavelength
> to frequency), then -- even though a plot of the two vectors looks
> exactly the same -- int_tabulated gives a different answer.
>
> IDL> x = reverse(x) & y = reverse(y)
> IDL> print,int_tabulated(x,y)
> 42.3619
>
> A closer look at the int_tabulated documentation tells us that if X if
> not monotonic increasing, then one needs to set the /SORT keyword.
> So now we follow the documentation:
>
> IDL> print,int_tabulated(x,y,/sort)
> 21.3333
>
> and get the right answer again. But as the documentation also "sort"
> of warns
> you, setting the /SORT keyword reorders the input X and Y variables.
> This is
> probably not a problem unless one has modified the Y vector on the fly
> by
> some factor f
>
> f = findgen(5) + 0.5
> print,int_tabulated(x, y*f,/sort)
>
> Upon output, the X vector has been sorted, but the Y vector has not,
> since it
> was passed as part of a temporary variable y*f. So one has lost
> the one to
> one correspondence between X and Y, screwing up all subsequent
> processing.
>
> --Wayne
|
|
|