Re: 2D Integration
- To: mathgroup at yoda.physics.unc.edu
- Subject: Re: 2D Integration
- From: keiper
- Date: Tue, 22 Jun 93 13:41:11 CDT
Robert Singleton wrote that multidimensional NIntegrate[ ] is inadequate:
> In my original post I specified that
> I was working with an InterpolatingFunction. In this case there is
> *NOTHING* that I have been able to do to make NIntegrate "decent" -
> except I now use Terry Robb's very nice routine ArrayIntegrate[]
> that integrates a list of data (or I use Maple).
ArrayIntegrate[ ] provides no estimate of the error and is non-adaptive.
You can get the same effect with NIntegrate[ ] by setting MinRecursion
and MaxRecursion:
In[15]:= f = Interpolation[Flatten[Table[{x, y, Sin[x^2 + x y + y^2]},
{x, 0., 2., .1}, {y, 0., 2., .1}], 1]]
Out[15]= InterpolatingFunction[{{0., 2.}, {0., 2.}}, <>]
In[16]:= g[n_] := Timing[NIntegrate[f[x, y], {x, 0, 2}, {y, 0, 2},
MinRecursion -> n, MaxRecursion -> n,
Compiled -> False, AccuracyGoal -> 0]]
In[17]:= g[0]
Out[17]= {0.566667 Second, 0.540633}
In[18]:= g[1]
Out[18]= {0.766667 Second, 0.545296}
In[19]:= g[2]
Out[19]= {3.48333 Second, 0.544633}
In[20]:= g[3]
Out[20]= {14.25 Second, 0.544647}
(Note that the timings increase quadratically with the recursion depth.)
Note that I have also set Compiled -> False. This saves a small amount of
time because the pseudo-complier is not (yet) able to deal with
lists and attempting to compile an InterpolatingFunction[ ] will just
waste time.
> > After all D[] knows how to handle InterpolatingFunction's, so why
> > not NIntegrate[] and Integrate[] too?
But ND[ ] does not do anything special with InterpolatingFunction[ ]s. The
point is well taken that Integrate[ ] should be able to undo anything that
D[ ] does, but I disagree that NIntegrate[ ] should do symbolic computation:
everytime we make NIntegrate[ ] "smarter" the overhead gets a little slower
and users wonder why something as trivial as NIntegrate[x^2, {x, 2, 5}]
takes so long. The "N" in NIntegrate means that the algorithm is numerical.
> This is VERY good suggestion! Why not make NIntegrate[] simply use
> the data points when it receives an InterpolatingFunction?!?
One very good reason is that if f is an InterpolatingFunction[ ] and you
want to integrate f[x^2] Sin[x] or even f[x^2] you have to do a LOT of
computation to either do a reparameterization of the integral or you
have to find the values of x that cause you to land on the data points.
This problem gets much worse in several dimensions where the grid of the
InterpolatingFunction[ ] may not even be aligned with the coordinate
axes.
There are three separate reasons why NIntegrate[ ] is not as fast as one
might like it to be:
1) The pseudo-compiler cannot deal with InterpolatingFunction[ ]s.
(This is being worked on.)
2) Evaluation of InterpolatingFunction[ ]s is not very fast.
3) The algorithms for NIntegrate[ ] assume that the integrand is
continuous in its first several derivatives. An
InterpolatingFunction[ ] is only continuous, but not smooth.
This problem would be somewhat less if the recursion in
NIntegrate[ ] caused the region of integration to be divided
up at the data points rather than always bisected. The
problem with doing this is the same as the problem of just
using the data points directly: except in trivial cases this
is much more difficult than it first appears.
Jerry B. Keiper
keiper at wri.com