# Integrating polynomial functions over polytopes with an add-on package

There is a Mathematica package to evaluate integrals over polytopes:

http://library.wolfram.com/infocenter/Books/3652/

In the documentation (Functions.nb file) I find:

ipoly::usage = "ipoly[f[x1, x2, I, xn], {x1, x2, I, xn}, {{a11, a12, I, a1n}, {a21, a22, I, a2n}, I, {aJ1, aJ2, I, aJn}}, {b1, b2, I, bJ}] is the n-dimensional integral of f[I] over a finite volume bounded by an n-dimensional convex polytope P. P is defined to be all points which satisfy the J inequalities: aj1 x1 + aj2 x2 + I + ajn xn <= bj, 1 <= j <= J. Input form is ipoly[f, x, {c1 <= c2, c3 <= c4, I}] where ci's are linear in x.";

I am trying a very simple example (similar to the one presented in AboutFunctions.nb: to integrate the function: f[x,y] = x + y over the polytope described by the set of inequalities {0 <= x, x <= 1 - y, -1 <= y, x + 2 y <= 2} .

ipoly[
x+y,
{x,y},
{0 <= x, x <= 1 - y, -1 <= y, x + 2 y <= 2}
]


ipoly[x + y, {x, y}, ..]

I am unable to understand how the output can be the result of an integral? It looks like it is just giving me back the input itself.

If someone knows how to use this function ipoly[...], please tell me.

• I haven't read/downloaded/tried the package, but I would guess that the first argument to ipoly needs to be a function that takes as many variables as you provide in the next argument, whereas you've input an expression. For your case above, try ipoly[Plus, {x,y}, ...]
– rm -rf
Mar 18, 2013 at 19:16

In Version 10 you can now create ImplicitRegions and Integrate over them:

region = ImplicitRegion[And @@ {0 <= x, x <= 1 - y, -1 <= y, x + 2 y <= 2}, {x, y}]


Now

Integrate[x + y, {x, y} ∈ region]


2/3

This package is from 1992. Much has happened in the land of Mathematica since then ;) You can achieve the same thing with the builtin Integrate / NIntegrate functions with the help of Boole for defining the region inside the polytope since Mathematica version 5.1. Let's take your example:

Integrate[
(x + y) Boole[
And @@ {0 <= x, x <= 1 - y, -1 <= y, x + 2 y <= 2}
],
{x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}
]
(* 2/3 *)


Also RegionPlot can be helpful to check the polytope region visually

RegionPlot[ And @@ {0 <= x, x <= 1 - y, -1 <= y, x + 2 y <= 2}, {x, -1, 3}, {y, -2, 2}]


• Thanks a lot for that quick help. It seems like its good for my purpose. However, I'll still wait for a while before am sure of it. Mar 18, 2013 at 19:19
• Sure, just experiment a bit with it and see if it fits your needs. Also have a look at the reference documentation i linked to. I has some good examples how to integrate over regions. Mar 18, 2013 at 19:23
• This method certainly works well for small number of inequality constraints. However, my application involves more constraints that I think mathematica can handle with the above implementation. Just to confirm that I had taken some time - I had 5 free variables and 10 liner inequalities. When I started the computation on nohup, it ran for a while and said: No More Memory Available, Mathematica Kernel has shut down. Solving the inequalities is the real memory intensive process. I was wondering if I could directly input the vertices of a polytope instead of inequalities, it might be faster. Mar 19, 2013 at 15:17
• Maybe if you update your question with an example of the kind of polytopes you are dealing with, more people could have a look and see what solutions there might be. Apr 1, 2013 at 22:07