6
$\begingroup$

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]] [[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
        ) // TraditionalForm
    ]

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.

$\endgroup$
3
  • 1
    $\begingroup$ a bit aside, but this is not an indefinite integral. $\endgroup$
    – george2079
    Jan 21, 2018 at 12:31
  • $\begingroup$ Regarding your final comment: it's not exactly $\LaTeX$, but MathJax. It is not as flexible as $\LaTeX$. $\endgroup$
    – anderstood
    Jan 21, 2018 at 15:49
  • $\begingroup$ 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[5]]] $\endgroup$
    – Carl Woll
    Jan 22, 2018 at 18:17

4 Answers 4

6
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ Very nice! But it would be helpful if you included the definition of deltaFunctionIntegral. $\endgroup$
    – Carl Woll
    Jan 22, 2018 at 18:58
  • 1
    $\begingroup$ Awesome! And it'll be better if you append some mathematical explanation(/deducing procedure) of the equation. Thx!! $\endgroup$
    – ooo
    Jan 23, 2018 at 4:06
  • $\begingroup$ Quite lucid !! Looking forward about some mathematical explanation too :))... $\endgroup$ Jan 23, 2018 at 11:11
  • 1
    $\begingroup$ 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.) $\endgroup$
    – ooo
    Jan 29, 2018 at 8:58
  • $\begingroup$ Sorry, I thought I'd done that, but apparently the edit didn't take. :( I'll try again later today. $\endgroup$
    – Pillsy
    Jan 29, 2018 at 16:29
6
$\begingroup$

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)$

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

(This answer is just a simplification of @Pillsy's answer)

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[99]] //AbsoluteTiming

{0.02938, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

By comparison, deltaFunctionIntegral takes about 1000 times longer:

deltaFunctionIntegral[100, y] //AbsoluteTiming

{31.0896, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

$\endgroup$
1
  • $\begingroup$ Wow!! You have given a very significant and essential improvement! $\endgroup$
    – ooo
    Jan 23, 2018 at 8:38
2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ I agree that sometimes ImplicitFunction is not the faster method. Even the generation of the vars is slower in my code. Voting up ! $\endgroup$ Jan 21, 2018 at 14:07
  • 1
    $\begingroup$ 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? $\endgroup$
    – ooo
    Jan 21, 2018 at 15:30
  • $\begingroup$ @ooo You might want to indicate the version and OS you're using. $\endgroup$
    – anderstood
    Jan 21, 2018 at 15:47
  • 1
    $\begingroup$ @JoséAntonioDíazNavas this has been observed before, for some reason using atomic symbols ( a1,a2,.. ) is faster than indexed symbols ( a[1],a[2] ) for these things. You can also do like var=Table[Unique[],{n}] $\endgroup$
    – george2079
    Jan 21, 2018 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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