       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:

Out= 5.55742
In:= data = Table[Exp[-x+2 y^2]//N,{x,-1,1,0.05},{y,0,1,0.05}];
Out= 5.57303
Out= -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

```

• Prev by Date: 0^t
• Next by Date: Computer Algebra Systems emphasizing Mathematica