Integrate returns "diverging integral" for obviously convergent expression

Question

Bug introduced in 13.0 or earlier. Fixed in 13.2.0 or earlier.

---

Integrate[Cos[phi]^2*Sin[phi]^n, {phi, 0, Pi}, Assumptions -> n > 0]
on Mathematica 13.0 returns the obviously mistaken error message (Integrate::idiv):

Curiously enough, just replacing $n$ by $n-2$,
Integrate[Cos[phi]^2*Sin[phi]^(n-2), {phi, 0, Pi}, Assumptions -> n > 2]
returns the correct answer $\frac{\sqrt{\pi }\, \Gamma \left(\frac{n-1}{2}\right)}{2 \Gamma \left(\frac{n}{2}+1\right)}.$

Any idea what is going on?

---

Thank you for the workarounds. Still, I would be interested to know why this fails, is there a lesson here, is it a more general difficulty? The integrand seems so simple and well-behaved, why does Mathematica think it "does not converge"?

Answer 1

I suggest a workaround

f[n_] = Integrate[Sin[phi]^n, {phi, 0, Pi}, Assumptions -> n > 0]
FullSimplify[f[n] - f[n + 2]]

$$\frac{\sqrt{\pi } \Gamma \left(\frac{n+1}{2}\right)}{2 \Gamma \left(\frac{n}{2}+2\right)}$$

Another way of calculating the integral

FullSimplify[
 Integrate[Cos[phi]^2*Sin[phi]^n, {phi, 0, Pi/4, Pi}, 
  Assumptions -> n > 0]]

suggests that in MA way of integrating the antiderivative has a singularity at $\phi=\pi/4$. Have no idea why.

Answer 2

Here's a way to get the answer by changing the integral by symmetry. Powers of sine and cosine that depend on parameters have always been tricky for Integrate. Not sure why the trick below works and the original fails.

Integrate[
 Cos[phi]^2*Sin[phi]^n + Sin[phi]^2*Cos[phi]^n, {phi, 0, Pi/2}, 
 Assumptions -> n > 0]
(*  (Sqrt[π] Gamma[(1 + n)/2])/(2 Gamma[2 + n/2])  *)

Workarounds from my comments:

Both

Integrate[Cos[phi]^2*Sin[phi]^n, {phi, 0, Pi}, 
  GenerateConditions -> False] // FullSimplify[#, n < -4] &

Normal@Integrate[#, {phi, 0, Pi}] & /@ 
  Expand[(1 - Sin[phi]^2) Sin[phi]^n] // FullSimplify

give

(*  (Sqrt[π] Gamma[(1 + n)/2])/(2 Gamma[2 + n/2])  *)

but not in a satisfying way.

And

Integrate[Cos[phi]^2*Sin[phi]^n // TrigToExp, {phi, 0, Pi}, Assumptions -> n > 0]

gives

(*  (π^(3/2) Sec[(n π)/2])/(2 Gamma[1/2 - n/2] Gamma[2 + n/2])  *)

which you have to use Limit to evaluate at odd integers. Or we can expand secant and gamma according to their infinite products to get the first result above:

(π^(3/2) Sec[(n π)/2])/(2 Gamma[1/2 - n/2] Gamma[2 + n/2]) /. 
  Sec[(n π)/2]/Gamma[1/2 - n/2] -> 1/(Pi ((π - n Pi)/(2 π))*
    Product[
     (1 - (π - n π)^2 /
      (4 k^2 π^2)) ((1 + 1/k)^(1/2 (-1 - n))/(1 + (-1 - n)/(2 k))),
     {k, Infinity}]) // FullSimplify

---

Update:

Integrate tries a substitution $\phi = \arctan(p/2)$, which fails for reasons I don't understand, since the transformed integral, while improper, converges:

Cos[phi]^2*Sin[phi]^n*Dt[phi, p] /. phi -> ArcTan[p/2] /. 
   t : _Sin | _Cos :> TrigExpand[t] /. p -> phi // Simplify
Integrate[2 %, {phi, 0, \[Infinity]}, Assumptions -> n > 0]
(*
  (8 (phi/Sqrt[4 + phi^2])^n)/(4 + phi^2)^2
  (Sqrt[\[Pi]] Gamma[(1 + n)/2])/(2 Gamma[2 + n/2])  
*)

My guess is that there's a simple(?) bug that is not particularly interesting, but maybe someone else will have more insight. Integrate can generate answers to equivalent integrals involving things like Hypergeometric2F1[2, 1/2 + n/2, 3/2, 1], which is undefined for n>0. It evaluates to ComplexInfinity if n is given a definite, nonnegative value. I do not know if that is connected to the idiv divergence error.

Answer 3

This seems to work:

res = Integrate[Cos[phi]^m*Sin[phi]^n, {phi, 0, \[Pi]}] /. m -> 2


N[Table[{NIntegrate[Cos[phi]^2*Sin[phi]^n, {phi, 0, Pi}], 
  res}, {n, 1, 6}]]

(* {{0.666667, 0.666667}, {0.392699, 0.392699}, {0.266667, 
 0.266667}, {0.19635, 0.19635}, {0.152381, 0.152381}, {0.122718, 0.122718}} *)

Answer 4

Two small steps are better than one big step. A simple workaround in 13.0.0 on Windows 10 is as follows.

Integrate[1/2*Cos[2*phi]*Sin[phi]^n, {phi, 0, Pi},  Assumptions -> n > 0] + 
Integrate[1/2*Sin[phi]^n, {phi, 0, Pi}, Assumptions -> n > 0]

(Sqrt[\[Pi]] Gamma[(1 + n)/2])/(2 Gamma[1 + n/2]) - ( n Sqrt[\[Pi]] Gamma[(1 + n)/2])/(4 Gamma[2 + n/2])

FullSimplify[%]

(Sqrt[\[Pi]] Gamma[(1 + n)/2])/(2 Gamma[2 + n/2])

Addition. Concerning possible explanation of the incorrect result, I share @MichaelE2's opinion since

Integrate[Cos[phi]^2*Sin[phi]^n, phi, Assumptions -> n > 0]

(1/(1 + n))Sqrt[Cos[phi]^2] Hypergeometric2F1[-(1/2), (1 + n)/2, (3 + n)/2, Sin[phi]^2] Sec[ phi] Sin[phi]^(1 + n)

Plot3D[%, {phi, 0, Pi}, {n, 0, 3}]