calculus - How to evaluate $int_{pi/2}^pi sqrt{1 - frac{1}{2}cos^2 x + sin x sin 2x} ;mathrm{d}x$

Here is an integral from a very old math journal:
(https://www.jstor.org/stable/1967417)

$$\int_{\pi/2}^\pi \sqrt{1 - \frac{1}{2}\cos^2 x + \sin x \sin 2x} \;\mathrm{d}x$$
It appears the journal never published a solution to this problem.

My attempts:

1. If we make the substitution $t = \cos x$, then the integral equals
   $$\int_{-1}^0 \sqrt{\frac{1 - \frac{t^2}{2} + 2t(1 - t^2)}{1 - t^2}} \;\mathrm{d}t$$
   However it appears impossible to integrate this function. Mathematica fails to evaluate both the indefinite and the definite integral.




2. If we integrate this function numerically, we get 0.827760002939144239418157727592. The Inverse Symbolic Calculator returns nothing for this number. WolframAlpha gives nothing interesting either.

Update: The integral can be expressed as a double sum. However it converges slowly. I am now looking for a faster way to compute its value.

Answer

It is actually very easy to express the integral as an infinite series. Note that
$$\sqrt{1 + x} = 1 + \sum_{k = 1}^\infty c_k x^k$$
$$c_k = (-1)^{k+1}\frac{(2k - 3)!!}{k! \cdot 2^k}$$
and
$$\sin x \sin 2x = 2\sin^2 x \cos x$$
Then

$$\begin{aligned}I &= \int_{\pi/2}^\pi\sqrt{1 + \sin x \sin 2x - \frac{1}{2}\cos^2 x} \;\mathrm{d}x\\ &= \int_{\pi/2}^\pi 1 + \sum_{k = 1}^\infty c_k \left(2\sin^2 x \cos x - \frac{1}{2} \cos^2 x\right)^k \;\mathrm{d}x\\ &= \frac{\pi}{2} + \sum_{k = 1}^\infty c_k \int_0^{\pi/2} \left(-2\cos^2 x \sin x - \frac{1}{2} \sin^2 x\right)^k \;\mathrm{d}x\\ &= \frac{\pi}{2} + \sum_{k = 1}^\infty (-1)^k c_k \int_0^{\pi/2} \sum_{m = 0}^k \binom{k}{m}\left(2\cos^2 x\sin x\right)^m\left(\frac{1}{2}\sin^2 x\right)^{k - m} \;\mathrm{d}x \end{aligned}$$
Now we have
$$\left(2\cos^2 x\sin x\right)^m\left(\frac{1}{2}\sin^2 x\right)^{k - m} = 2^{2m - k} \sin^{2k - m} x \cos^{2m} x$$
Table of Integrals, Series, and Products (8th edition, p. 397, formula 3.621-5) tells me that
$$\int_0^{\pi/2} \sin^m x \cos^n x \;\mathrm{d}x = \frac{1}{2}\mathrm{B}\left(\frac{m + 1}{2},\frac{n + 1}{2}\right)$$

where $\mathrm{B}(a,b)$ is the complete Beta function. Therefore
$$I = \frac{\pi}{2} + \sum_{k = 1}^\infty (-1)^k c_k \sum_{m = 0}^k \binom{k}{m} 2^{2m - k - 1} \mathrm{B}\left(\frac{2k-m+1}{2},\frac{2m+1}{2}\right)$$
The numerical results seem to agree, but the series converges very slowly. The first 500 terms sum up to 0.82776000293935595965, which still has an error of $10^{-12}$.

Because this is a double sum, we need $O(n^2)$ time to calculate the first $n$ terms. I have not figured out a faster way to calculate.

Update: This problem oddly reminds me of a binomial sum problem on ProjectEuler (Problem 567). Note that the Beta function satisfies
$$\mathrm{B}(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a + b)}$$
For multiples of $1/2$, we can compute the Gamma function explicitly
$$\Gamma\left(\frac{n}{2}\right) = \left\{\begin{aligned} \left(\frac{n}{2} - 1\right)! & \quad n \text{ even}\\ \sqrt{\pi}\frac{(n - 2)!!}{2^\frac{n - 1}{2}} & \quad n \text{ odd}\end{aligned}\right.$$

Let us rewrite the integral as
$$I = \frac{\pi}{2} - \sum_{k = 1}^\infty \sum_{m = 0}^k d_{k,m}$$
$$d_{k,m} = \frac{(2k - 3)!!}{k! \cdot 2^k} \binom{k}{m} 2^{2m - k - 1} \mathrm{B}\left(\frac{2k-m+1}{2},\frac{2m+1}{2}\right)$$
The calculation of the Beta function term will depend on the parity of $m$. It is easy to verify that
$$\frac{d_{k,2t + 2}}{d_{k,2t}} = \frac{4(k - 2t - 1)(k - 2t)(4t + 1)(4t + 3)}{(2k - 2t - 1)(k + t + 1)(t + 1)(2t + 1)}$$
$$\frac{d_{k,2t + 3}}{d_{k,2t + 1}} = \frac{4(k - 2t - 2)(k - 2t - 1)(4t + 3)(4t + 5)}{(k - t - 1)(2k + 2t + 3)(t + 1)(2t + 3)}$$
$$\frac{d_{k + 1,2t}}{d_{k,2t}} = \frac{(2k - 2t + 1)(2k - 1)}{8(k - 2t + 1)(k + t + 1)}$$

$$\frac{d_{k + 1,2t + 1}}{d_{k,2t + 1}} = \frac{(2k - 1)(k - t)}{2(k - 2t)(2k + 2t + 3)}$$
From here I think it is possible to transform the double sum into a recurrence relation, using hypergeometric summation techniques. I will post results if I have time to do it.

Update: I tried Zeilberger's algorithm. It works but the results look very messy, and is probably not a good idea to get the numerical value of this series. Since Zeilberger's algorithm is a late 20th century invention, I do not think the original problem poser has a better way to calculate the integral, other than adding the series by hand. See Zeilberger's book (https://www.math.upenn.edu/~wilf/AeqB.html) for a description of the algorithm if you want to repeat the attempt.

Update again: I am very impressed by the numerical capabilities of Mathematica. It turns out that we can simply compute
$$I = E + O$$
$$E = \sum_{t = 0}^\infty d_{2t,2t} \cdot {}_2F_1\left(t + \frac{1}{2},2t - \frac{1}{2};3t + 1;\frac{1}{2}\right)$$
$$O = \sum_{t = 0}^\infty d_{2t + 1,2t + 1} \cdot {}_2F_1\left(t + 1,2t + \frac{1}{2};3t + \frac{5}{2};\frac{1}{2}\right)$$
Using Mathematica I summed up the first 450 terms in each series. The result is 0.827760002939144239418158294687, which has an error of $10^{-23}$. I think this is the best we can extract out of this problem. I will stop working on it.