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