4
$\begingroup$

hopefully someone can help me with this problem I am having. I am trying to calculate the integral

Integrate[(((a - 2) (2 x^2 - 1/2) + 2)) (((T x (x - 1)) + M)^(-a/2)), {x, 0, 1},  
GenerateConditions -> False]

Which should equal

Integrate[(((a - 2) (2 x^2 - 1/2) + 2)) (((T x (x - 1)))^(-a/2)), {x, 0, 1},
 GenerateConditions -> False]

when I afterwards set M to zero. However, what I get is that the first integral equals

(2^(-3 + a) (10 - 3 a) Sqrt[\[Pi]] (-T)^(-a/2) Gamma[1 - a/2])/Gamma[5/2 -
 a/2]

Which is correct, yet in the limit M to 0, the first integral returns "ComplexInfinity". I have tried using variations of assumptions as well as no assumptions but I always encounter the same problem. I am relatively new to Mathematica, so I am not sure whether this is a bug or I am doing something wrong. Any help that people can offer is greatly appreciated.

$\endgroup$
  • $\begingroup$ No matter what the sign of T, one of the last two arguments of AppellF1 is infinite at M == 0. Instead, try expanding the first solution about M == 0. $\endgroup$ – bbgodfrey Apr 13 '17 at 19:59
  • $\begingroup$ The two results are not equal, as can be demonstrated by evaluating them at {T -> 1, a -> 1, M -> 1/8}. However, there is no reason to suppose that they should be, because they may be valid only for different conditions on the parameters. $\endgroup$ – bbgodfrey Apr 13 '17 at 20:19
  • $\begingroup$ bbgodfrey, the solutions should be the same when M is specified to be zero surely, as then the two integrands become equal. If the two functions are being defined on different domains is there any way to specify the domains to be the same? I have tried adding the assumptions T, M, a all real as well as M,a >=0, but this does not change that the two integrals give different results. Also expanding around M =0 is not ideal as i need to take the limit a ->0 in the end, and that could introduce new problems $\endgroup$ – Ragemeistro Apr 13 '17 at 21:17
  • 1
    $\begingroup$ You can use NIntegrate in specific cases (i.e., for specific numerical values of the parameters) to verify which formulas are correct. $\endgroup$ – Michael E2 Apr 13 '17 at 23:15
  • 2
    $\begingroup$ I think you misunderstood. I was not suggesting you use NIntegrate in place of Integrate. The main weakness of my suggestion is that it can be used only as a spot-check to show a formula is wrong at a point value of {T, M, a}, but that would still be informative. $\endgroup$ – Michael E2 Apr 14 '17 at 10:41
4
$\begingroup$

In situations like this, it is prudent to obtain from Integrate the conditions under which a solution it provides is "guaranteed" valid. This is accomplished by setting the option GenerateConditions to True instead of False. Then

int = (((a - 2) (2 x^2 - 1/2) + 2)) (((T x (x - 1)) + M)^(-a/2));
Integrate[int, {x, 0, 1}, GenerateConditions -> True]

produces a ConditionalExpression a bit long (LeafCount of 696) to reproduce here. However,

FullSimplify[% /. M -> 0]

produces a less lengthy and more informative expression

(* ConditionalExpression[1/6 (0^(-a/2) 2^(2 + a) (-2 + a) 
   AppellF1[3, a/2, a/2, 4, (2 T)/(T - Sqrt[T^2]), (2 T)/(T + Sqrt[T^2])] + 
   (3 (-6 + a) T (Beta[(T - Sqrt[T^2])/(2 T), 1 - a/2, 1 - a/2] - 
                  Beta[(T + Sqrt[T^2])/(2 T), 1 - a/2, 1 - a/2]))/SqrtT^2]), 
   Re[T/Sqrt[T^2]] > 1 || Re[T/Sqrt[T^2]] < -1 || T/Sqrt[T^2] ∉ Reals] *)

Because this condition is not satisfied by any complex number, the result with M == 0 is not valid, not does Integrate claim that it is. So, the behavior described in the question is not a bug. Rather, Integrate is unable to obtain a solution for general M that is also valid for M == 0. This conclusion is corroborated by the fact that

Integrate[int, {x, 0, 1}, GenerateConditions -> True, Assumptions -> M >= 0]

returns unevaluated. Additional insight can be obtained by considering a specific case,

f[M0_?NumericQ] := NIntegrate[int /. {a -> 1, T -> 1, M -> M0}, {x, 0, 1}, 
    MaxRecursion -> 20]
Plot[{Re@f[M], Im@f[M]}, {M, -1, 1}, PlotRange -> {-7, 10}]

enter image description here

Evidently, the integral is not differentiable with respect to M at M == 0, although it is continuous there. Integrate gives the same curve with

Integrate{int /. {a -> 1, T -> 1}, {x, 0, 1}, GenerateConditions -> True, 
    Assumptions -> M >= 0]
(* ConditionalExpression[-Sqrt[M] + (7/2 + 2 M) ArcCoth[2 Sqrt[M]], M > 1/4] *)
Plot[Evaluate@ReIm @(% // First), {M, -1, 1}, Exclusions -> None, PlotRange -> {-7, 10}]

even for M < 1/4. (Whether an overly restrictive condition constitutes a bug is, perhaps, a matter of opinion.) On the other had, omitting the Assumption yields

Integrate[int /. {a -> 1, T -> 1}, {x, 0, 1}, GenerateConditions -> True]
(* ConditionalExpression[1/4 (-4 Sqrt[M] + (7 + 4 M) Log[-1 - 2 Sqrt[M]] - 
   (7 + 4 M) Log[1 - 2 Sqrt[M]]), Im[Sqrt[1 - 4 M]] != 0 || Re[Sqrt[1 - 4 M]] > 1] *)

For real M, Integrate asserts that the answer is valid when M < 0, and indeed this is true.

Plot[Evaluate@ReIm @(%13 // First), {M, -1, 1}, Exclusions -> None, PlotRange -> {-7, 10}]

enter image description here

Results with Assumptions on Constants

The OP recently provided constraints on the three constants appearing in the integrand, M >= 0, T < 0, and 0 <= a < 2. Additionally, it is convenient to represent M as m T, with m <= 0. Then, the integrand becomes

int = T^(-a/2) ((a - 2) (2 x^2 - 1/2) + 2) (x (x - 1) + m)^(-a/2)

with T now entering only through a multiplicative factor. Applying the constraints then yields

Integrate[int, {x, 0, 1}, Assumptions -> 0 <= a < 2 && m < 0, GenerateConditions -> True];
sn = FullSimplify[%, 0 <= a < 2 && m < 0]
(* 1/6 T^(-a/2) (4 (-2 + a) E^(-(1/2) I a π) (-(1/m))^(a/2)
   AppellF1[3, a/2, a/2, 4, 2/(1 + Sqrt[1 - 4 m]), -(2/(-1 + Sqrt[1 - 4 m]))] - 
   3 I (-6 + a) m^(1 - a/2) (-((I m)/Sqrt[1 - 4 m]))^(1/2 (-2 + a)) (-1 + 4 m)^(-a/4)
   (Beta[1/2 - 1/(2 Sqrt[1 - 4 m]), 1 - a/2, 1 - a/2] - 
   Beta[1/2 (1 + 1/Sqrt[1 - 4 m]), 1 - a/2, 1 - a/2])) *)

Results can be obtained for positive m as well, although the results are a bit long to be reproduced here.

s4 = Integrate[int, {x, 0, 1}, Assumptions -> 0 <= a < 2 && 0 < m < 1/4, GenerateConditions -> True];
sp = Integrate[int, {x, 0, 1}, Assumptions -> 0 <= a < 2 && m > 1/4, GenerateConditions -> True];

Note that Integrate returned no additional constraints for any of these three results. Consequently, the solution for all real m is

s = Piecewise[{{sn, m < 0}, {sp, m > 1/4}}, s4]

with, as a sample result,

Plot[s /. {a -> 1, T -> -1}], {m, -1, 1}, Exclusions -> None, PlotRange -> {-10, Automatic}]

enter image description here

The corresponding numerical integration yields an identical plot.

f[m0_?NumericQ] := NIntegrate[int /. {a -> 1, T -> -1, m -> m0}, {x, 0, 1}, 
    MaxRecursion -> 20]
Plot[{Re@f[m], Im@f[m]}, {m, -1, 1}, PlotRange -> {-10, Automatic}]

s can, of course, be plotted for other values of a, although the evaluation of s is quite slow for some a with 0 < m < 1/4. Since the question requested only m < 0, the slow speed for 0 < m < 1/4 does not really matter.

$\endgroup$
  • $\begingroup$ Thanks for this, it is really detailed which is great! The only thing that confuses me is that there does exist a solution when 'M=0', so does this mean that Integrate is returning a solution just valid on a subdomain of the whole M-domain of the solution? That is to say the process that Integrate uses does not seem to be very general, but is there any way to make Integrate return an answer which is also valid when M=0? the restriction 'M>=0' does not seem to change the output of integrate. $\endgroup$ – Ragemeistro Apr 23 '17 at 11:25
  • 1
    $\begingroup$ @Ragemeistro Precisely so. Integrate seems unable to obtain a solution valid for all real M. The option GenerateConditions -> True instructs Integrate to provide the range of parameters for which it knows its solution to be correct. However, I have found a more general solution including M in the vicinity of zero, but I need an answer to my comment of a few minutes ago before proceeding. $\endgroup$ – bbgodfrey Apr 23 '17 at 17:15

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.