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