MathGroup Archive 1993

[Date Index] [Thread Index] [Author Index]

Search the Archive

2D Integration

  • To: mathgroup at
  • Subject: 2D Integration
  • From: Robert Singleton <bobs at>
  • 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.


bob singleton 


bobs at

  • Prev by Date: 0^t
  • Next by Date: Computer Algebra Systems emphasizing Mathematica
  • Previous by thread: 0^t
  • Next by thread: Re: 2D Integration