# Examples

➤ 1.

$$I=\iiint\limits_D \dfrac{{\rm d}x\ \!{\rm d}y\ \!{\rm d}z}{(x+y+z)^p}~,$$

where $$D=\{~(x,y,z)~\colon~x>0~,y>0~,z>0~,x+y+z<1~\}~.$$

MMA gives: $~I=-\dfrac1{2(p-3)}~.$

➤ 2.

$$I=\idotsint\limits_{~~V} \dfrac{{\rm d}x_1\ \!{\rm d}x_2\ \!\cdots\ \!{\rm d}x_n}{(x_1+x_2+...+x_n)^p}~,$$

where $$V=\{~(x_1,x_2~,\dots,x_n)~\colon~x_1>0~,x_2>0~~,\dots,x_n>0~,x_1+x_2+...+x_n<1~\}~.$$

The integral converges when $p<n$.

Verified the case of $n=1～5$, conjecture by induction that $$I=\frac{1}{(n-1)!~(n-p)}~.$$

# Conclusion

Special thanks to @Pillsy's formula and @CarlWoll's simplification!

### In:

n = 20;
Integrate[y^(n-1-p), {y,0,1}, x∈Simplex[n-1]] [] // AbsoluteTiming


### Out:

$\left\{0.363959,\frac1{121645100408832000 (20-p)}\right\}$

NOTE: The following content is no longer worth browsing.

# Code

### In:

int[var_] :=
Block[
var,
fun = Power[Plus @@ var, p]~Power~-1;

con0 = <|1 -> (StringRiffle[# > 0 & /@ var, "&&"] // ToExpression),
2 -> Plus @@ var < 1|>;

con = StringRiffle[Values@con0, "&&"] // ToExpression;

Column[{fun, con}] // TraditionalForm // Print;

(result = Integrate[fun, var ∈ ImplicitRegion[con, Evaluate@var],
Assumptions -> p < Length@var] // Together
]

int@{x, y, z}
int@{x, y, z, u, v, w} // Timing


### Out:

$(x+y+z)^{-p} \\ x>0\land y>0\land z>0\land x+y+z<1 \\ -\dfrac1{2(p-3)}$

$(u+v+w+x+y+z)^{-p} \\ x>0\land y>0\land z>0\land u>0\land v>0\land w>0\land u+v+w+x+y+z<1 \\ \left\{29.6719,-\dfrac1{120(p-6)}\right\}$

# Timing

$\begin{array}{rr} n & t(s) \\ \hline 1 & 0.328 \\ 2 & 0.656 \\ 3 & 1.391 \\ 4 & 3.297 \\ 5 & 16.828 \\ 6 & 29.641 \\ 7 & 69.766 \end{array}$

Tests on:
Mathematica 11.1.1.0 (←2017-05-16),
Windows 64-bit 10.0.10586 (←2015-11-12),
AMD A8-5600K (←2012-09-26)

# Aiming

Code of fast and efficient algorithm.

Wish: regardless of how large is $n$, the computation time is always less than 10s.

I would be thankful, if you could improve or rewrite any part of the code.

I would be admired, if you could write a library/macro that can output very fast and can apply to many situations, whereas even do not depend on the original Integrate function.

Question END.

The following has nothing to do with the post topic, just recording my experiences of fixing this —

Your post appears to contain code that is not properly formatted as code.

In my case, all the problems stem from LaTeX MathJax (comment from @anderstood).
① Custom commands are not allowed to submit (though they are rendered normally).
② Need to delete after \\ in some places.

• a bit aside, but this is not an indefinite integral. Jan 21 '18 at 12:31
• Regarding your final comment: it's not exactly $\LaTeX$, but MathJax. It is not as flexible as $\LaTeX$. Jan 21 '18 at 15:49
• You can use Simplex for your region specification, but it won't be much (if any) faster, e.g. With[{v=Table[Unique[x],{5}]},Integrate[1/Total[v]^p, v\[Element]Simplex]] Jan 22 '18 at 18:17

I was able to speed things up incredibly by observing that the integrand only depends on the sum $x_1 + x_2 + \ldots + x_n$, and we can take advantage of this by introducing a gratuitous integration over a new variable and throwing in a Dirac $\delta$ "function", which satisfies the following relation:

$$\int_a^b dx\,f(x)\,\delta(x-c) = \begin{cases} f(c) & a < c < b \\ 0 & c < a \text{ or } b < c \end{cases}$$

With a little coaxing, you can reproduce this in Mathematica, using DiracDelta:

Integrate[f[x]*DiracDelta[x - c], {x, a, b},
Assumptions -> {a < b, #}] & /@ {
a < c < b,
c < a || b < a
}
(* {f[c], 0} *)


Now, use a specific application of the defining relation. For $0 < x_1 + \ldots + x_n < 1$, we have

$$\frac{1}{\left(x_1 + \ldots + x_n\right)^{\,p}} = \int_0^1 dy \frac{\delta\left(y - (x_1 + \ldots + x_n\right)}{y^p}$$

We can insert this into the original expression, and swap around the order of the integration.

$$\idotsint\limits_{V} \frac{dx_1\cdots\,dx_n}{(x_1 + \ldots + x_n)^p } = \int_0^1 \frac{dy}{y^p} \int_0^\infty dx_1 \dots \int_0^\infty dx_n \,\,\, \delta\left(y - (x_1 + \ldots + x_n\right)$$

Now we can do the integral over the $x_i$s separately

deltaFunctionIntegral[n_Integer?Positive, y_] :=
With[{vars = Array[\[FormalX], n]},
Apply[
Integrate[DiracDelta[y - Total[vars]], ##,
Assumptions -> Positive[y]] &,
Map[{#, 0, Infinity} &, vars]]];


Mathematica will blaze through it:

Table[deltaFunctionIntegral[n, y] // AbsoluteTiming, {n, 1, 10}]
(*
{{0.024917, 1},
{0.033459, y},
{0.042226, y^2/2},
{0.045033, y^3/6},
{0.050282, y^4/24},
{0.05671, y^5/120},
{0.063595, y^6/720},
{0.073157, y^7/5040},
{0.080864, y^8/40320},
{0.086789, y^9/362880}}
*)


Now we can do the integral with this factor and things remain speedy:

Table[
With[{factor = deltaFunctionIntegral[n, y]},
Integrate[factor/y^p, {y, 0, 1}]] // AbsoluteTiming
{n, 1, 10}]
(*
{{0.288039, ConditionalExpression[1/(1 - p), Re[p] < 1]},
{0.300244, ConditionalExpression[1/(2 - p), Re[p] < 2]},
{0.314402, ConditionalExpression[1/(2 (3 - p)), Re[p] < 3]},
{0.32506, ConditionalExpression[1/(6 (4 - p)), Re[p] < 4]},
{0.318678, ConditionalExpression[1/(24 (5 - p)), Re[p] < 5]},
{0.339789, ConditionalExpression[1/(120 (6 - p)), Re[p] < 6]},
{0.346237, ConditionalExpression[1/(720 (7 - p)), Re[p] < 7]},
{0.35325, ConditionalExpression[1/(5040 (8 - p)), Re[p] < 8]},
{0.365873, ConditionalExpression[1/(40320 (9 - p)), Re[p] < 9]},
{0.367103, ConditionalExpression[1/(362880 (10 - p)), Re[p] < 10]}}
*)


I know how to show that deltaFunctionIntegral[n, y] gives $y^{n-1}/(n-1)!$ outside of Mathematica for general $n$, but I don't know how to coax the result out of Mathematica yet.

• Very nice! But it would be helpful if you included the definition of deltaFunctionIntegral. Jan 22 '18 at 18:58
• Awesome! And it'll be better if you append some mathematical explanation(/deducing procedure) of the equation. Thx!!
– ooo
Jan 23 '18 at 4:06
• Quite lucid !! Looking forward about some mathematical explanation too :))... Jan 23 '18 at 11:11
• Oh~ Can you please make some mathematical explanation on that integral of Dirac function? (I can't understand even if I try to give examples.)
– ooo
Jan 29 '18 at 8:58
• Sorry, I thought I'd done that, but apparently the edit didn't take. :( I'll try again later today. Jan 29 '18 at 16:29

This is what I get in my 7-years old Intel I3 in my macOS 10.13.2 and MMA ver 11.2:

Table[AbsoluteTiming[vars = Array[x, n];
Integrate[(Apply[Plus, vars])^-p, vars \[Element]
ImplicitRegion[And @@ (# > 0 & /@ vars) && Apply[Plus, vars] < 1,Evaluate@vars]]], {n, 7}]


$\left( \begin{array}{cc} 0.488669 & \text{ConditionalExpression}\left[\frac{1}{1-p},\Re(p)<1\right] \\ 0.978526 & \text{ConditionalExpression}\left[\frac{\frac{1}{2-p}-1}{p-1},\Re(p)<2\right] \\ 2.07194 & \text{ConditionalExpression}\left[\frac{\frac{1}{3-p}-\frac{p}{2}}{p^2-3 p+2},\Re(p)<3\right] \\ 4.77045 & \text{ConditionalExpression}\left[\frac{p^2-2 p-\frac{6}{4-p}+3}{-6 p^3+36 p^2-66 p+36},\Re(p)<4\right] \\ 20.0045 & \text{ConditionalExpression}\left[-\frac{1}{24 (p-5)},\Re(p)<5\right] \\ 33.4052 & \text{ConditionalExpression}\left[-\frac{1}{120 (p-6)},\Re(p)<6\right] \\ 84.0427 & \text{ConditionalExpression}\left[\frac{\frac{720}{7-p}-p ((p-7) p ((p-7) p+28)+252)}{720 (p-6) (p-5) (p-4) (p-3) (p-2) (p-1)},\Re(p)<7\right] \\ \end{array} \right)$

• Thanks! Your rewriting is much better than mine!
– ooo
Jan 22 '18 at 8:35
• However, slower than that from @george2079 :((... Anyway, always happy to contribute and help Jan 22 '18 at 8:36
• I wondering my tests result (comment to @george2079). So, can you show a single test on $n=7$ compare yours with his?
– ooo
Jan 22 '18 at 8:42
• Ummm..., mine 80.6 sec, and that of @george2079 101.5 in a clear kernel... Jan 22 '18 at 8:52
• ( ＾◡＾)っ The "self-optimization" of MMA has not been defeated yet~
– ooo
Jan 22 '18 at 9:10

You could make the change of variables $x_i \to y x_i$ to convert the integral to $\underset{x\in \Delta^{n-1}}{\int } y^{n-1} dx$ (here $\Delta^{n-1}$ is the usual nomenclature for Simplex[n-1]) which evaluates very quickly. For example:

Integrate[y^99, x ∈ Simplex] //AbsoluteTiming


{0.02938, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

By comparison, deltaFunctionIntegral takes about 1000 times longer:

deltaFunctionIntegral[100, y] //AbsoluteTiming


{31.0896, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

• Wow!! You have given a very significant and essential improvement!
– ooo
Jan 23 '18 at 8:38

Old school, directly construct the integration limits instead of using a region proves to be just slightly faster.

n = 5
var = Symbol["a" <> ToString[#]] & /@ Range[n]
fun = 1/(Total@var)^p
Integrate[fun, Evaluate@Sequence @@ MapThread[{#1, 0, 1 - #2} &,
{var, Accumulate[Prepend[var[[;; -2]], 0]]}],
Assumptions -> p < Length@var] // Timing


{9.90606, 1/(120 - 24 p)}

n=6


{29.5466, 1/(720 - 120 p)}

It doesn't seem to make much difference in timing if you leave out the assumption on p and let it generate a conditional.

• I agree that sometimes ImplicitFunction is not the faster method. Even the generation of the vars is slower in my code. Voting up ! Jan 21 '18 at 14:07
• Running on my AMD A8-5600K, as $n$ increasing, it is very strange that yours use more time. AbsoluteTiming's stable values — Yours: {n, t} = {5, 10.5} | {6, 36.6} | {7, 160}. Mine: {n, t} = {5, 21.6} | {6, 33.8} | {7, 103}. — Can anyone explain what on earth happened inside the MMA compiler for this example?
– ooo
Jan 21 '18 at 15:30
• @ooo You might want to indicate the version and OS you're using. Jan 21 '18 at 15:47
• @JoséAntonioDíazNavas this has been observed before, for some reason using atomic symbols ( a1,a2,.. ) is faster than indexed symbols ( a,a ) for these things. You can also do like var=Table[Unique[],{n}] Jan 21 '18 at 16:06