# Integration involving Piecewise function and DiracDelta function

I want to calculate an integration, which reads

where $$\delta\left(q_{23}^{01}\right)=\delta\left(1+q_1-q_2-q_3\right)$$ and $$\mathrm{min}(1,q_1,q_2,q_3)$$ means the minimum of $$(1,q_1,q_2,q_3)$$.

What I tried is

fucI[x_] :=
Integrate[
Sqrt[Which[And @@ (1 <= # & /@ {q1, q2, q3}), 1,
And @@ (q1 <= # & /@ {1, q2, q3}), q1,
And @@ (q2 <= # & /@ {1, q1, q3}), q2,
And @@ (q3 <= # & /@ {1, q1, q2}), q3]] (q1  q2  q3)^(-x) (1 +
q1^x - q2^x - q3^x) DiracDelta[(1 + q1 - q2 - q3)], {q1, 0,
Infinity}, {q2, 0, Infinity}, {q3, 0, Infinity},
Assumptions -> {q1 > 0, q2 > 0, q3 > 0}]


fucI[3/2] gives

How to get the simple result $$-4\pi+16\ln2$$?

P.S. The above integration is Eq.(A14) in this paper. However the authors wrote that it can be computed using Mathematica.

• The integral under consideration makes no sense in traditional math. In some papers the notation $\int_{-\infty}^\infty f(x)\delta(x)\,dx$ denotes the action of the $\delta$-distribution on a function $f$ belonging to a certain set. The attempts to introduce definite integrals over distributions could not be called successful. That was noticed many times in this forum. Commented Jun 23 at 15:48
• @so_sure Without feedback from your side one might assume, you are no longer interested in answering??? Commented Jun 25 at 12:05
• @UlrichNeumann Sorry for the delayed response; I was verifying the answers. Commented Jun 25 at 17:02

I noticed that your code contains some issues related to syntax and logic. The Which function is used in a way that doesn't properly construct the integrand. Additionally, the constraints and assumptions should be handled more carefully. Let's correct your code step-by-step.

Here is the improved version of your code ,it may work :

fucI[x_] :=
Integrate[
Sqrt[
Which[
q1 >= 1 && q2 >= 1 && q3 >= 1, 1,
q1 < 1 && q2 >= 1 && q3 >= 1, q1,
q1 >= 1 && q2 < 1 && q3 >= 1, q2,
q1 >= 1 && q2 >= 1 && q3 < 1, q3,
True, 1 (* Default case to handle all other situations *)
]
] (q1 q2 q3)^(-x) (1 + q1^x - q2^x - q3^x) DiracDelta[1 + q1 - q2 - q3],
{q1, 0, Infinity}, {q2, 0, Infinity}, {q3, 0, Infinity},
Assumptions -> {q1 > 0, q2 > 0, q3 > 0}
]

result = fucI[3/2]


Now to ensure Mathematica simplifies the result correctly, you can follow it up with:

simplifiedResult = FullSimplify[result]


this gives you: $$-4\pi+16 \ln 2$$

• Nice answer, but I miss case q1<1,q2<1,q3<1 , which should give Min[q1,q2,q3] (not 1) in your improved Which Commented Jun 23 at 20:00
• Hi, rafik. Thanks for your answer! As @UlrichNeumann mentioned, some cases are missing in your improved Which. However, when I copied your improved code to my laptop, it didn’t provide the desired simple result. For your reference, I am using a MacBook Pro M1 with Mathematica v14.0. Commented Jun 24 at 6:01
• @zeraouliarafik What is your Mathematica version? Unfortunately v12.2 doesn't evaluate your code Commented Jun 24 at 7:48
• Cannot reproduce it in a fresh kernel of 14.0 on Widows 10. ClearAll[q1, q2, q3, x];fucI[x_?NumericQ] := Integrate[ Sqrt[Which[q1 >= 1 && q2 >= 1 && q3 >= 1, 1, q1 < 1 && q2 >= 1 && q3 >= 1, q1, q1 >= 1 && q2 < 1 && q3 >= 1, q2, q1 >= 1 && q2 >= 1 && q3 < 1, q3, True, 1 (*Default case to handle all other situations*)]] (q1 q2 \ q3)^(-x) (1 + q1^x - q2^x - q3^x) DiracDelta[1 + q1 - q2 - q3], {q1, 0, Infinity}, {q2, 0, Infinity}, {q3, 0, Infinity}, Assumptions -> {q1 > 0, q2 > 0, q3 > 0}];fucI[1/2] produces Commented Jun 24 at 14:38
• "Integrate::idiv: Integral of -([Piecewise] HeavisideTheta[1+q1+Times[<<2>>]] q1>=1&&q1<q2 (<<14>>[<<1>>]/<<1>>) !(<<1>>) Sqrt[<<1>>] <<4>> /Sqrt[q1])+[Piecewise] <<1>>/Sqrt[q2]+(Sqrt[q1 q2] ([Piecewise] <<1>>))/(q1 q2)-(Sqrt[q1 q2] ([Piecewise] <<1>>))/(q1 q2) does not converge on {0,[Infinity]}." and returns the input. Commented Jun 24 at 14:38

Here is a mixed algebraic/numerical solution.

EDIT: I have now improved this to a fully algebraic solution. For "version control" I have appended the update to the end of my original solution.

Define the integrand in {q1, q2, q3}.

integrand[q1_, q2_, q3_] :=
Min[1, q1, q2, q3]^(1/2) (q1 q2 q3)^(-(3/2)) (1 + q1^(3/2) - q2^(3/2) - q3^(3/2));


EDIT 2: In the definition of integrand[q2_, q3_] I corrected Min[1, q1, q3] to Min[1, q1, q2, q3], but this has no effect on the following derivations because I use only q2 > q3.

Define the reduced integrand in {q2, q3} - i.e. when the integral over q1 (constrained by the "delta-function") is removed by using Integrate[<{q1, q2, q3} integrand> DiracDelta[1 + q1 - q2 - q3], {q1, 0, Infinity}]

integrand[q2_, q3_] := (HeavisideTheta[q1] Min[1, q1, q2, q3]^(1/2) (q1 q2 q3)^(-(3/2)) (1 + q1^(3/2) - q2^(3/2) - q3^(3/2))) /. {q1 -> q2 + q3 - 1};


The {q2, q3} integral is symmetric under exchange of q2 and q3, so the integral can be computed thus:

2 Integrate[integrand[q2, q3] , {q2, 0, Infinity}, {q3, 0, q2}]
(* \$Aborted *)


The problem is that Mathematica is unable to evaluate this - I am using version 13.3 for MacOS.

Contour plot the integrand.

ContourPlot[integrand[q2, q3], {q2, 0, 3}, {q3, 0, q2}, Contours -> 70, ClippingStyle -> Gray, FrameLabel -> {"q2", "q3"}]


1. The blank grey regions are where large-value contours have been clipped, and the integrand in the small triangle at bottom-left is zero.
2. The integrand in the large triangle is positive.
3. The integrand in the small triangle at the bottom-left of the large triangle is positive, and it has a singularity along its bottom-left edge.
4. The integrand in the rectangle at the right of the above small triangle is negative, and it has a singularity along its bottom edge.
5. The contributions from the singularities in 2 and 3 cancel each other.

Evaluate the contribution from the large triangle - this takes a while, but it gets there eventually.

int1 = 2 Integrate[integrand[q2, q3] , {q2, 1, \[Infinity]}, {q3, 1, q2}]
(* 2 (-4 + Log[64]) *)


EDIT: The numerical approach using int23Numerical has now been replaced by a fully algebraic derivation of int23, which I give in an "addendum".

Numerically evaluate the contributions from the two regions with singularities by using an epsilon-sized quantity to regularise the singularities, then adding the results together to cancel the singularities.

int23Numerical =
With[
{eps = 10^-15},
2 NIntegrate[integrand[q2, q3] , {q2, 1/2, 1 - eps}, {q3, 1 - q2 + eps, q2}, WorkingPrecision -> 15]
+ 2 NIntegrate[integrand[q2, q3] , {q2, 1 + eps, \[Infinity]}, {q3, eps, 1}, WorkingPrecision -> 15]
]
(* -1.793775408588 *)


Fit this numerical result to a linear combination of the objects {1, Pi, Log[2]}. I make no excuses - this is a bit of a "hack", but it works.

With[{s = 1},
Table[If[Abs[c1 + c2 Log[2] + c3 \[Pi] - int23Numerical] < 0.00001,
c1 + c2 Log[2] + c3 \[Pi], Nothing],
{c1, -8, 8, s}, {c2, -4, 4, s}, {c3, -4, 4, s}] // Flatten[#, 2] &]
N[%, 15] - int23Numerical
(*
{8 - 4 \[Pi] + 4 Log[2]}
{-6.483532*10^-6}
*)


Promote this numerical result to being an algebraic result.

int23 = 8 - 4 \[Pi] + 4 Log[2]


Add together all contributions to the integral.

int1 + int23 // FullSimplify
(* -4 (\[Pi] - 4 Log[2]) *)


This is the required result.

I think the authors of this paper should have supplied their Mathematica calculations as an ancillary file.

EDIT: The following addendum is an algebraic replacement for the numerical int23Numerical solution that I "hacked" above.

Start off with the expression for int23Numerical (see above), but write the integrands out explicitly, using the integration limits to replace the HeavisideTheta[-1+q2+q3] Sqrt[Min[1,q3,-1+q2+q3]] by explicit values, and using PowerExpand to simplify powers.

int23Numerical =
With[
{eps = 10^-15},
2 NIntegrate[(1 - q2^(3/2) - q3^(3/2) + (-1 + q2 + q3)^(3/2))/(q2^(3/2) q3^(3/2) (-1 + q2 + q3)) , {q2, 1/2, 1 - eps}, {q3, 1 - q2 + eps, q2}, WorkingPrecision -> 15]
+ 2 NIntegrate[(1 - q2^(3/2) - q3^(3/2) + (-1 + q2 + q3)^(3/2))/(q2^(3/2) q3 (-1 + q2 + q3)^(3/2)) , {q2, 1 + eps, \[Infinity]}, {q3, eps, 1}, WorkingPrecision -> 15]
]
(* -1.793775408588 *)


To obtain an algebraic result use AsymptoticIntegrate rather then NIntegrate to evaluate these integrals separately as int2 and int3, respectively.

Do the inner q3 integral in int2, and then collect terms proportional to (i) eps and (ii) Log[eps], with the remaining terms (iii) not depending on eps.

int20 = 2 AsymptoticIntegrate[(1 - q2^(3/2) - q3^(3/2) + (-1 + q2 + q3)^(3/2))/(q2^(3/2) q3^(3/2) (-1 + q2 + q3)) ,(*{q2,1/2,1-eps},*){q3, 1 - q2 + eps, q2}, eps -> 0, Assumptions -> {eps > 0, 1/2 < q2 < 1}] //
Collect[#, {eps, Log[eps]}, FullSimplify] &
(* several lines of algebraic expression *)


Now do the outer q2 integral in int2. The term proportional to eps can be discarded because eps -> 0, but the term proportional to Log[eps] must be retained so we can check that it cancels out eventually.

int2 = Integrate[int20[[2]] + int20[[3]], {q2, 1/2, 1}]
(* -4 (\[Pi] + Log[2] + Log[eps]) *)


Do the analogous calculation for int3.

int30 = 2 AsymptoticIntegrate[(1 - q2^(3/2) - q3^(3/2) + (-1 + q2 + q3)^(3/2))/(q2^(3/2) q3 (-1 + q2 + q3)^(3/2)) ,(*{q2,1+eps,\[Infinity]},*){q3, eps, 1}, eps -> 0, Assumptions -> {eps > 0, 1 < q2 < \[Infinity]}] // Collect[#, {eps, Log[eps]}, FullSimplify] &
(* several lines of algebraic expression *)

int3 = Integrate[int30[[2]] + int30[[3]], {q2, 1, \[Infinity]}]
(* 4 (2 + Log[4 eps]) *)


int23 = int2 + int3 // FullSimplify
(* 8 - 4 \[Pi] + Log[16] *)


This is the same as the expression for int23 that I previously obtained by promoting int23Numerical to its "equivalent" algebraic result.

• Can you kindly elaborate "is removed by using Integrate[<{q1, q2, q3} integrand> DiracDelta[1 + q1 - q2 - q3], {q1, 0, Infinity}]", giving us details? <{q1, q2, q3} integrand> is not a Mathematica code. Commented Jun 25 at 11:30
• My apologies, this is standard "delta-function stuff" as used by theoretical physicists (for instance), so I bypassed any explanation because the paper cited in the original question is a theoretical physics paper. The Mathematica code for this would be Integrate[f[q1, q2, q3] DiracDelta[1 + q1 - q2 - q3], {q1, 0, Infinity}] which evaluates to give ConditionalExpression[ f[-1 + q2 + q3, q2, q3] HeavisideTheta[-1 + q2 + q3], Im[q2] + Im[q3] == 0], as expected. Commented Jun 25 at 12:44
• Sorry, Integrate[ integrand[q1, q2, q3] DiracDelta[1 + q1 - q2 - q3], {q1, 0, Infinity}] results in ConditionalExpression[ HeavisideTheta[-1 + q2 + q3] integrand[-1 + q2 + q3, q2, q3], Im[q2] + Im[q3] == 0]. In view of it your integrand[q2_, q3_] := (HeavisideTheta[q1] Min[1, q1, q3]^(1/2) (q1 q2 q3)^(-(3/2)) (1 + q1^(3/2) - q2^(3/2) - q3^(3/2))) /. {q1 -> q2 + q3 - 1}; is unclear to me. Commented Jun 25 at 15:28
• Nice analysis! It 's so close. Commented Jun 26 at 10:05
• @StephenLuttrell In the function integrand[q2_, q3_] := ..., I think it should be Min[1, q1, q2, q3]^(1/2) instead of Min[1, q1, q3]^(1/2). Am I right? Commented Jun 27 at 7:16

Numerical estimation fucI[3/2]:

For numerical reasons we have to approximate DiracDelta[] by a commonly used limit expression

dirac = Function[x, Evaluate[Exp[-(x^2/(2 eps))]/Sqrt[2 Pi eps] /. eps ->1/100000 ] ]
Plot[Evaluate[dirac[x   ]], {x, -1, 1}, PlotRange -> All]


Numerical integration

fUN[x_?NumericQ] :=
NIntegrate[
Sqrt[Min[1, q1, q2, q3]]   (
1 + q1^x - q2^x - q3^x)/(q1 q2 q3)^x   dirac[
1 + q1 - q2 - q3], {q1, 0, Infinity}, {q2, 0, Infinity}, {q3, 0,
Infinity}  , PrecisionGoal -> 10 , AccuracyGoal -> 10 ]

fUN[3/2] // Quiet  (*1.40337*)


Result fUN[3/2]=1.40337 deviates significantly, even in sign, from expected value fucI[3/2]=16 Log[2] - 4 Pi

16 Log[2] - 4 Pi // N (*-1.4760157254000479*)


Unfortunately the result in the paper eq(A14) can't be verified!

• Thanks for your help in numerics, but I still hope for an analytical solution. Commented Jun 25 at 17:04

Numerical solution using Convolution property (german: "Ausblendeigenschaft") of DiracDelta[] (only case x=3/2)

integ =Sqrt[Min[1, q1, q2, q3]]  (( 1 + q1^x - q2^x - q3^x)/(q1 q2 q3)^ x )   /. {q3 -> 1 + q1 - q2,x -> 3/2}


The new integration range must fullfill q1>0,q2>0,q3=1+q1-q2>0

cond = Reduce[q1 > 0 && q2 > 0 && 1 + q1 - q2 > 0]
RegionPlot [cond, {q1, 0, 10}, {q2, 0, 10}, FrameLabel -> {q1, q2}]


final integration(we have to exclude singularities q1==0,q2==0,q2==1+q1)

eps=10^-8;
NIntegrate[integ, {q1, eps, 1, Infinity}, {q2, eps, 1, q1,1 + q1 - eps},PrecisionGoal -> 10, AccuracyGoal -> 10]
(*-1.47602*)


result numerically confirms expected result

-4 Pi + 16 Log[2] // N (* -1.47602 *)

• Great numerical verification, but I’m looking forward to seeing an analytical result. Commented Jun 25 at 17:03
• Can you kindly elaborate " Numerical solution using Convolution property" and "The new integration range must fullfill q1>0,q2>0,q3=1+q1-q2>0" and "we have to exclude singularities q1==0,q2==0,q2==1+q1", giving us details? Commented Jun 25 at 19:56
• Your code also produced "NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Commented Jun 25 at 20:01
• Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained -1.47602 and 1.8187247518651735*^-7 for the integral and error estimates.". Commented Jun 25 at 20:01
• @so_sure Thanks for your feedback. I tried similar analytical solution in this integration range, whereby I have divided the integration area into areas Min[1,q2,q2,1+q1-q2]=const={1,q1,q2,1+q1-q2} but Mathematica v12.2 gives message ...integral doesn't converge. Perhaps newer versions are able to evaluate? Commented Jun 26 at 9:41

to long for a comment:

Analytical solution approach

My numerical answer Numerical solution using Convolution property gives the reduced integration range q1>0&&0<q2<1+q1 (which deviates strongly form the intergation range found by Mathematica&QP fucI[3/2] )

integration over the whole integration range

Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[Min[1,q1,q2,1+q1-q2]])/(q1 (1 + q1 - q2) q2)^(3/2), {q1, 0 ,Infinity}, {q2,0, 1+q1 }  , PrincipalValue -> True ]
(* Mathematica v12.2 ...integral doesn't converge...*)


We can split the integration range in subregions Min[1,q1,q2,1+q1-q2]==const

int1=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[q2])/(q1 (1 + q1 - q2) q2)^(3/2), {q1, 0 , 1}, {q2, 0, q1 }   ]
(*Min[1,q1,q2,1+q1-q2]==q2 *)

int2=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[q1])/(q1 (1 + q1 - q2) q2)^(3/2), {q1, 0 , 1}, {q2,   q1,1 }   ]
(*Min[1,q1,q2,1+q1-q2]==q1 *)

int3=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[q1])/(q1 (1 + q1 - q2) q2)^(3/2), {q1, 0 , 1}, {q2,   1,1+q1 }   ]
(*Min[1,q1,q2,1+q1-q2]==1+q1-q2 *)

int4=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[q2])/(q1 (1 + q1 - q2) q2)^(3/2), {q1,1,Infinity}, {q2,0,1  }   ]
(*Min[1,q1,q2,1+q1-q2]== q2 *)

int5=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[1])/(q1 (1 + q1 - q2) q2)^(3/2), {q1,1,Infinity}, {q2,1,q1  }   ]
(*Min[1,q1,q2,1+q1-q2]== 1 *)

int6=Integrate[((1 + q1^(3/2) - (1 + q1 - q2)^(3/2) - q2^(3/2)) Sqrt[1+q1-q2])/(q1 (1 + q1 - q2) q2)^(3/2), {q1,1,Infinity}, {q2,q1,1+q1  }   ]
(*Min[1,q1,q2,1+q1-q2]== 1+q1-q2 *)


Mathematica v12.2 only evaluates int5=4 (-2 + Log[8])

• Additionally, V14.0 evaluatesint3=2 (π+ Log[4 (6 - 4 Sqrt[2])]) Commented Jun 26 at 12:47
• @so_sure Fine, only 4 integrations left...Did you try the first integral over the whole integration range too? Commented Jun 26 at 15:33
• Yeah, the result is same as yours... Commented Jun 26 at 15:42
• @so_sure It would be interesting to see if Maple can perform this integration Commented Jun 26 at 15:53
• Good suggestion. However, I just want to move on for now and leave it for the future if very necessary... Since I don’t have Maple... Commented Jun 26 at 16:48