2D Integration

*To*: mathgroup at thnext.mit.edu*Subject*: 2D Integration*From*: Robert Singleton <bobs at thnext.mit.edu>*Date*: Tue, 1 Jun 93 14:26:17 -0400

A while back I posted a question on 2D integration - i.e. is there a good way to do this in Mma? It would be nice if WRI would make a decent 2D or 3D integration routine (hint hint) -- the built in Mma NIntegrate is much to slow for multiple integration involving interpolating functions. Until then I had hoped to use InterCall and pass an interpolating function into a fortran or C program that will do the integration. I'm impressed on the whole with InterCall but unfortunately, for this problem it doesn't quite do the trick (passing the data to the fortran program is too slow -- InterCall cannot yet compile an interpolating function and so it has to continually ask Mma to evaluate the function.). I do have an acceptable solution to the problem and people might be interested in it. Terry Robb wrote the following nifty Mma routines that integrate lists of data. They are very fast - the time is spent in generating the data itself. Here is a nice routine for 2D and 3D integration: ListIntegrate1D[data_List, h_] := h*(Apply[Plus,data] - (First[data] + Last[data])/2); ListIntegrate2D[data_List, {dx_, dy_}] := ListIntegrate1D[Map[ListIntegrate1D[#,dy]&, data],dx]; ListIntegrate3D[data_List, {dx_, dy_, dz_}] := ListIntegrate1D[Map[ListIntegrate2D[#,dy,dz]&, data],dx]; Here is a test run of ListIntegrate2D[] -- I haven't tested the 3D routine: In[10]:= answer1 = NIntegrate[Exp[-x+2 y^2],{x,-1,1},{y,0,1}] Out[10]= 5.55742 In[11]:= data = Table[Exp[-x+2 y^2]//N,{x,-1,1,0.05},{y,0,1,0.05}]; In[12]:= answer2 = ListIntegrate2D[data, {0.05,0.05}] Out[12]= 5.57303 In[13]:= (answer1-answer2)/answer2 Out[13]= -0.00280218 This is about 0.3% accurate -- very good. Here is an even better solution (again from Terry Robb): ListIntegrate[data_List, dx_] := (Apply[Plus,data] - (First[data] + Last[data])/2) * dx; ListIntegrate[data_List, dx_, dydz__] := ListIntegrate[Map[ListIntegrate[#,dydz]&, data],dx]; and here is a compact way of doing cubic interpolation: ArrayIntegrate[data_List, dx_, dydz__] := ArrayIntegrate[Map[ArrayIntegrate[#,dydz]&, data], dx]; ArrayIntegrate[data_List, dx_] := (* trapezium rule *) (Apply[Plus, data] - (First[data] + Last[data])/2)*dx; ArrayIntegrate[data_List, dx_] := (* Simpson's rule *) (Apply[Plus, Partition[data,2]].{2,4} - First[data] + Last[data] )/3*dx /; OddQ[Length[data]]; These routines do n-dim integration but I haven't tested them yet. Hope this helps some people. Regards, bob singleton MIT bobs at thnext.mit.edu