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$
  • 1
    $\begingroup$ a bit aside, but this is not an indefinite integral. $\endgroup$ – george2079 Jan 21 '18 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 '18 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 '18 at 18:17
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$
  • 1
    $\begingroup$ Very nice! But it would be helpful if you included the definition of deltaFunctionIntegral. $\endgroup$ – Carl Woll Jan 22 '18 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 '18 at 4:06
  • $\begingroup$ Quite lucid !! Looking forward about some mathematical explanation too :))... $\endgroup$ – José Antonio Díaz Navas Jan 23 '18 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 '18 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 '18 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$
  • $\begingroup$ Thanks! Your rewriting is much better than mine! $\endgroup$ – ooo Jan 22 '18 at 8:35
  • $\begingroup$ However, slower than that from @george2079 :((... Anyway, always happy to contribute and help $\endgroup$ – José Antonio Díaz Navas Jan 22 '18 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 '18 at 8:42
  • $\begingroup$ Ummm..., mine 80.6 sec, and that of @george2079 101.5 in a clear kernel... $\endgroup$ – José Antonio Díaz Navas Jan 22 '18 at 8:52
  • $\begingroup$ ( ^◡^)っ The "self-optimization" of MMA has not been defeated yet~ $\endgroup$ – ooo Jan 22 '18 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$
  • $\begingroup$ Wow!! You have given a very significant and essential improvement! $\endgroup$ – ooo Jan 23 '18 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$
  • $\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$ – José Antonio Díaz Navas Jan 21 '18 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 '18 at 15:30
  • $\begingroup$ @ooo You might want to indicate the version and OS you're using. $\endgroup$ – anderstood Jan 21 '18 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 '18 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.