Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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 '13 at 19:16
add comment

1 Answer

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}]

Polytope RegionPlot

share|improve this answer
    
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. –  Pavithran Iyer Mar 18 '13 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. –  Thies Heidecke Mar 18 '13 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. –  Pavithran Iyer Mar 19 '13 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. –  Thies Heidecke Apr 1 '13 at 22:07
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.