"Integral Milking" I begin this post with a plea: please don't be too harsh with this post for being off topic or vague. It's a question about something I find myself doing as a mathematician, and wonder whether others do it as well. It is a soft question about recreational mathematics - in reality, I'm shooting for more of a conversation.
I know that a lot of users on this site (e.g. Cleo, Jack D'Aurizio, and so on) are really good at figuring out crafty ways of solving recreational definite integrals, like
$$\int_{\pi/2}^{\pi} \frac{x\sin(x)}{5-4\cos(x)} \, dx$$
or
$$\int_0^\infty \bigg(\frac{x-1}{\ln^2(x)}-\frac{1}{\ln(x)}\bigg)\frac{dx}{x^2+1}$$
When questions like this pop up on MSE, the OP provides an integral to evaluate, and the answerers can evaluate it using awesome tricks including (but certainly not limited to):


*

*Clever substitution

*Exploitation of symmetry in the integrand

*Integration by parts

*Expanding the integrand as a series

*Differentiating a well-know integral-defined function, like the Gamma or Beta functions

*Taking Laplace and Inverse Laplace transforms


But when I play around with integrals on my own, I don't always have a particular problem to work on. Instead, I start with a known integral, like
$$\int_0^\pi \cos(mx)\cos(nx) \, dx=\frac{\pi}{2}\delta_{mn},\space\space \forall m,n\in
\mathbb Z^+$$
and "milk" it, for lack of a better word, to see how many other obscure, rare, or aesthetically pleasing integrals I can derive from it using some of the above techniques. For example, using the above integral, one might divide both sides by $m$, getting
$$\int_0^\pi \frac{\cos(mx)}{m}\cos(nx) \, dx=\frac{\pi}{2m}\delta_{mn},\space\space \forall m,n,k\in
\mathbb Z^+$$
Then, summing both sides from $m=1$ to $\infty$, and exploiting a well-known Fourier Series, obtain
$$\int_0^\pi \cos(nx)\ln(2-2\cos(x)) \, dx=-\frac{\pi}{n},\space\space \forall n\in
\mathbb Z^+$$
or, after a bit of algebra, the aesthetically pleasing result
$$\int_0^{\pi/2} \cos(2nx)\ln(\sin(x)) \, dx=-\frac{\pi}{4n},\space\space \forall n\in
\mathbb Z^+$$
After pulling a trick like this, I look through all of my notebooks and integral tables for other known integrals on which I can get away with the same trick, just to see what integrals I can "milk" out of them in the same way. This is just an example - even using the same starting integral, countless others can be obtained by using other Fourier Series, Power Series, integral identities, etc. For example, some integrals derived from the very same starting integral include
$$\int_0^\pi \frac{\cos(nx)}{q-\cos(x)} \, dx=\frac{\pi(q-\sqrt{q^2-1})^{n+1}}{1-q^2+q\sqrt{q^2-1}}$$
$$\int_0^\pi \frac{dx}{(1+a^2-2a\cos(x))(1+b^2-2b\cos(mx))}=\frac{\pi(1+a^m b)}{(1-a^2)(1-b^2)(1-a^m b)}$$
and the astounding identity
$$\int_0^{\pi/2}\ln{\lvert\sin(mx)\rvert}\cdot \ln{\lvert\sin(nx)\rvert}\, dx=\frac{\pi^3}{24}\frac{\gcd^2(m,n)}{mn}+\frac{\pi\ln^2(2)}{2}$$
Everyone seems to be curious about the proof of this last identity. A proof can be found in my answer here.
I just pick a starting integral, and using every technique I know as many times as possible, try to come up with the most exotic integrals as I can, rather than picking a specific integral and trying to solve it.
Of course, integrals generated this way would be poor (or at least extremely difficult) candidates for contest problems or puzzles to evaluate given the integral, since they are derived "backwards," and determining the derivation given the integral is likely much harder than pursuing the vague goal of a "nice-looking integral" with no objective objective (ha ha).
QUESTION: Do you (residents of MSE who regularly answer/pose recreational definite integral questions) do this same activity, in which you try to generate, rather than solve, cool integrals? If so, what are some  integrals you have come up with in this way? What strategies do you use? Does anyone care to opine on the value (or perhaps lack of value) of seeking integrals in this way?
Cheers!
 A: Yes, especially definite integrals. To begin, I often use the function $\ln$ and its counterpart $\exp$ to start 'building' the integrand. Then I plug it in Wolfram to check that the indefinite integral doesn't have a closed form (to make sure the integral is challenging enough :)
Seeing as there is no closed form, I use substitution/by parts to create another integral. The final step is to combine the integrals at the LHS to give a constant on the right.
Some examples:
$$\small\int_0^\infty\ln\left(\frac{e^x+1}{e^x-1}\right)\,dx=2\int_0^\infty\frac{xe^x}{e^{2x}-1}\,dx=2\int_1^\infty\frac{\ln x}{x^2-1}\,dx=\frac12\int_0^\infty\frac{\ln(x+1)}{x\sqrt{x+1}}\,dx=\frac12\int_0^\infty\frac{xe^{x/2}}{e^x-1}\,dx$$
and starting with the integral $\int_0^1\frac{\ln x}{e^x}\,dx$, we find that
$$\int_0^1\frac{(x^2-3x+1)\ln x}{e^x}\,dx=-\frac1e$$ 
The same method can be used to show that $$\int_0^\infty\frac{(x^2-3x+1)\ln x}{e^x}\,dx=0.$$
A: What a delightful thread!
I once spent many weeks playing with collections of line sources, free-surface pressure patches and other hydrodynamic singularities to estimate the pressures they induced on the sea-bed in finite-depth water.
Asked to give an impromptu talk to some naval engineers who wanted to know what I was working on at the time, I explained that E. O. Tuck had been using a very simple series to approximate a multiplicative term required in some integrals, and copied some notes I had to a blackboard:
$$2 \sum^{\infty}_{n=1} \frac{1}{(2n-1)^3} = \sum^{\infty}_{n=1} \frac{L_{n}}{n^{2}} = \frac{7}{4}\zeta(3)$$
where
$$L_n = 1 + \frac{1}{3} + \frac{1}{5} + \dots + \frac{1}{2n-1}.$$
The first sum is well known (23.3.20, A&S 1972, p. 807), and as there were
nods of approval from the engies I thought I'd just bore them with how the 2nd series arose, so I mumbled on for a few minutes more about other work we were doing.
I remember being very impressed with their easy familiarity with the zeta 
function. Older now, I'm not quite so sure that I read that audience response correctly. It was a Friday afternoon, last day of term, so I suspect that some had given the free departmental sherry a bit of a nudge after our long lunch.
A: Unsure if this is worthy of an answer, but one particular trick I find fascinating is coordinate changes that leave the result of an integration untouched. 
For example, there's a theorem with a name I can't remember right now (EDIT: It's called Glasser's master theorem, as Chappers pointed out below) that establishes equivalence of integrals of real functions over the entire real line: 
$$\int_{-\infty}^{\infty}f(x)dx = \int_{-\infty}^{\infty}f\left(|\alpha|x - \sum_{i=1}^{n}\frac{|\gamma_i|}{x - \beta_{i}}\right)dx$$
for arbitrary constants $\alpha$, $\beta_i$, $\gamma_i$. 
The reason why this is great for "milking" integrals is that you can keep changing the coordinates over and over until you get a monstrosity that has a simple result.
Let's try the simplest example I can think of, $\int_{-\infty}^{\infty}\frac{1}{x^2 + a}dx$ with real positive $a$. Then, by applying the coordinate change over and over using $\alpha = 1$, $\gamma_{i} = \gamma =1$ and $\beta_i = \beta = 0$:
$$\frac{\pi}{\sqrt{a}}=\int_{-\infty}^{\infty}\frac{1}{x^2 + a}dx = \int_{-\infty}^{\infty}\frac{x^2}{x^4 + (a+2)x^2 + 1}dx = \int_{-\infty}^{\infty}\frac{(x^2 (x^2 + 1)^2)}{a x^6 + 2 a x^4 + a x^2 + x^8 + 6 x^6 + 11 x^4 + 6 x^2 + 1} dx= \quad...$$
I'm certainly not suggesting this is a difficult integral, but you can see how it can get very hairy if I had nonzero $\beta_i$'s or more than one $\gamma_i$!
Once you have enough of these types of transforms under your belt, you can apply them to your heart's content in whatever form you like knowing that the result remains unchanged.
Hope this helps your quest to find more wonky integrals!
A: I'm not a very experienced integrator, but I've had some fun that I'd like to share. The below are results I've reached myself, but haven't been able to prove with complete confidence.
Let $\alpha_i \in\Bbb R, \forall i\in I=\{1,2,...,n\}$ when $n\in\Bbb N$, and $\forall i,j\in I, \alpha_i=\alpha_j \iff i=j$. Then:$$\int\prod_{i\in I}(\frac{1}{x-\alpha_i})dx=C+\sum_{i\in I}\frac{\ln |x-\alpha_i|}{\prod_{i≠j\in I}(\alpha_i-\alpha_j)}$$ Via partial fraction decomposition.
A: Mathematicians milk all sorts of things just for the fun of it!

*

*One Jon concocted a sequence of integrals using Fourier analysis that at first seem to evaluate to $\pi/2$ and then break down at the seventh term. He reported it as a bug to a computer algebra system vendor, who spent 3 days before figuring out he had been had.


*That whole MO thread above is about coming up with counter-examples to computer algebra systems such as difficult-to-see convergence. I also do that do Wolfram Alpha all the time; it will forever remain possible to easily trick it to give the Wrong Answer.


*The 1968 Putnam competition featured an obviously positive integral that evaluates to $\frac{22}{7}-\pi$. Lucas concocted another similar integral for $\frac{355}{113}-\pi$.


*You could ask the creator of $\int_{-\infty}^{+\infty} {dx \over 1 + \left(x + \tan x\right)^2} = \pi$ how he/she came up with this integral. It does not look like a naturally occurring integral to me.
A: We've seen Glasser's Master Theorem from ghostinthefigures well Ramanujan had one too:
Ramanujan's Master Theorem: Ramanujan was especially fond of this technique of his. G. H. Hardy stated that he "was particularly fond of them, and used them as one of his commonest tools (to milk his integrals)." Note that his procedure is strictly formal with results being valid only under severe conditions. 
Let $F(x)$ be some complex-valued function in some neighbourhood of $x=0$ with expansion
$$F(x)=\sum_{k=0}^{\infty}\frac{\lambda(k)(-x)^k}{k!}$$
where $\lambda$ is some analytic singular valued function. Then
$$
\int_{0}^{\infty}x^{n-1}F(x)\,dx=\Gamma(n)\lambda(-n)\tag{1}
$$
Alongside his Master Theorem $(1)$ Ramanujan also gave the equivalent identity:
$$\int _{0}^{\infty }x^{s-1}(\varphi(0)-x\varphi(1)+x^{2}\varphi(2)-\cdots)\,dx=\frac {\pi }{\sin(\pi s)}\varphi(-s)\tag{2}$$
which gets converted to the above form $(1)$ after substituting $\varphi(n)=\frac {\lambda(n)}{\Gamma (1+n)}$ and using the functional equation for the gamma function.
Ramanujan's proof starts with Euler's gamma function integral
$$\int_{0}^{\infty}e^{-mx}x^{n-1}\,dx=m^{-n}\Gamma(n),\quad m,\,n>0$$ 
Let $m=r^k$, $r$ a constant $>0$, and multiply both sides by $f^{(k)}(a)h^k/k!$, where $f$ is some function to be defined later with $a$, $h$, constants. Now sum on $k$ thus:
$$
\sum_{k=0}^{\infty}\frac{f^{(k)}(a)h^k}{k!} \int_{0}^{\infty}e^{-r^{k} x}x^{n-1}\,dx
=\Gamma(n)\sum_{k=0}^{\infty}\frac{f^{(k)}(a)(hr^{-n})^k}{k!}
$$
Now he expands $e^{-r^{k} x}$ in its Maclaurin series, inverts the order of summation and integration to get
\begin{align*}
\sum_{k=0}^{\infty}\frac{f^{(k)}(a)h^k}{k!} & \int_{0}^{\infty} x^{n-1}\sum_{j=0}^{\infty}\frac{r^{jk}(-x)^j}{j!}\,dx
= \int_{0}^{\infty} x^{n-1}\sum_{j=0}^{\infty}\frac{(-x)^j}{j!}\sum_{k=0}^{\infty}\frac{f^{(k)}(a)h^k r^{jk}}{k!}\,dx\\
&=\int_{0}^{\infty}x^{n-1}\sum_{j=0}^{\infty}\frac{f(a+hr^{j})(-x)^j}{j!}\,dx=\Gamma(n)f(a+hr^{-n})\tag{3}
\end{align*}
For $m$ real, let $f(hr^{m}+a)=\lambda(m)$, then $(3)$ can be written as $(1)$.
As an example of $(1)$ consider the binomial expansion
$$(1+x)^{-a}=\sum_{k=0}^{\infty}\binom{k+a-1}{k}(-x)^k =\sum_{k=0}^{\infty}\frac{\Gamma(k+a)}{\Gamma(a)}\frac{(-x)^k}{k!}$$
Let $\lambda(k)=\Gamma(k+a)/\Gamma(a)$. Then Ramanujan’s Master Theorem gives
$$\int_{0}^{\infty} \frac{x^{n-1}}{(1+x)^a}\,dx=\frac{\Gamma(n)\Gamma(a-n)}{\Gamma(a)}=B(n,a-n)$$
where $B$ is the beta integral.
As an example of $(2)$ consider the infinite product definition of the Gamma function:
$$\Gamma (x)=\frac {e^{-\gamma x}}{x}\prod _{n=1}^{\infty}\left(1+\frac {x}{n}\right)^{-1}e^{x/n}$$
is equivalent to the expansion
$$\log \Gamma (1+x)=-\gamma x+\sum _{k=2}^{\infty }\frac {\zeta (k)}{k}(-x)^{k}$$
where $\zeta (k)$ is the Riemann zeta function. This rearranges to
$$\frac{\gamma x+\log \Gamma (1+x)}{x^2}=\sum _{k=0}^{\infty }\frac {\zeta (2+k)}{2+k}(-x)^{k}$$
Let $\varphi(k)=\zeta (2+k)/(2+k)$. Then applying Ramanujan’s Master Theorem we have:
$$\int _{0}^{\infty }x^{s-1}\frac {\gamma x+\log \Gamma (1+x)}{x^{2}}\,dx=\frac {\pi }{\sin(\pi s)}\frac {\zeta (2-s)}{2-s}$$
valid for $0<\Re(s)<1$.
(See Ramanujan's Notebooks Vol I, Ch.4 for some of Ramanujan's results, and Hardy's Ramanujan Ch.XI, Definite Integrals, for thorough proofs.)
A: A correctly re-done version of my deleted answer.
For $a_k\ne 0$, $k=1,2,...,n$, we have
$$\Pi_n:=u\left(\small{\prod_{k=1}^{n}}a_k\right)\prod_{k=1}^{n}u(a_k)=u(1)+\sum_{k=1}^{n}\sum_{\sigma\subset S_n\\ |\sigma|=k}u\left(\small{\prod_{j\in\sigma}}a_j^2\right),\tag{1}$$ where $S_n=\{1,2,...,n\}$
(from here), with $u(x)=x+1/x$. Letting $a_k=\exp(i\theta_k)$, we have $u(a_k)=2\cos\theta_k$ and
$$u\left(\small{\prod_{j\in\sigma}}a_j^2\right)=u\left(\small{\prod_{j\in\sigma}}e^{2i\theta_j}\right)=u\left(e^{2i\small{\sum_{j\in\sigma}}\theta_j}\right)=2\cos\left(2\small{\sum_{j\in\sigma}}\theta_j\right).$$
Thus
$$\Pi_n=2^{n+1}\cos\left(\small{\sum_{k=1}^n\theta_k}\right)\prod_{k=1}^{n}\cos\theta_k=2+2\sum_{k=1}^{n}\sum_{\sigma\subset S_n\\ |\sigma|=k}\cos\left(2\small{\sum_{j\in\sigma}}\theta_j\right).$$
Setting $\theta_k=x$,
$$\Pi_n=2^n\cos(nx)\cos(x)^n=1+\sum_{k=1}^{n}\sum_{\sigma\subset S_n\\ |\sigma|=k}\cos(2|\sigma|x)=1+\sum_{k=1}^{n}{n\choose k}\cos(2kx),$$
because $${n\choose k}=\left|\left\{\sigma\subset S_n:|\sigma|=k\right\}\right|.$$
Thus
$$\int_0^{a}\cos(nx)\cos(x)^ndx=\frac1{2^n}\left(a+\sum_{k=1}^{n}{n\choose k}\frac{\sin(2ka)}{2k}\right),$$
which, I don't doubt, can be used to develop more interesting identities.

An identity developed from the above:
$$\int_0^1t\sin\left(n(n+1)t\right)\cos(nt)^{n-1}dt=\frac1{2^n n^2}\sum_{k=1}^{n}{n\choose k}\left(\frac{\sin(2nk)}{2nk}-\cos(2nk)\right),$$
where $n\in\Bbb Z_{\ge1}$.
A: Yes, definitely. For example, I found that
$$ m\int_0^{\infty} y^{\alpha} e^{-y}(1-e^{-y})^{m-1} \, dy = \Gamma(\alpha+1) \sum_{k \geq 1} (-1)^{k-1} \binom{m}{k} \frac{1}{k^{\alpha}} $$
(and related results for particular values of $\alpha$) while mucking about with some integrals. Months later, I was reading a paper about a particular regularisation scheme (loop regularisation) useful in particle physics, and was rather surprised to recognise the sum on the right! I was then able to use the integral to prove that such sums have a particular asymptotic that was required for the theory to actually work as intended, which the original author had verified numerically but not proved. The resulting paper's on arXiv here.
Never let it be said that mucking about with integrals is a pointless pursuit!
A: I'm not sure if this is what you are looking for, but I found the following (elliptic?) integral on AoPS:
$$\int \frac{x}{\sqrt{x^4+4x^3-6x^2+4x+1}}\, \mathrm{d}x$$
It was likely found by studying the Risch algorithm, or at least the Wikipedia page. Here it can be seen the general form of this integral cannot be solved in elementary terms, but by carefully choosing the right polynomials, one can construct a function whose derivative is of that form.
As shown here, the solution to the integral is the monstrous
$$ -\frac{1}{6}\ln\left[(x^{4}+10x^{3}+30x^{2}+22x-11)\sqrt{x^{4}+4x^{3}-6x^{2}+4x+1}- \\(x^{6}+12x^{5}+45x^{4}+44x^{3}-33x^{2}+43)\right]+C$$
This stems from the integral given in the Wikipedia page,
$$\int \frac{x}{\sqrt{x^4 + 10 x^2 - 96 x - 71}} \mathrm{d}x$$
which evaluates to 
$$ - \frac{1}{8}\ln \,\Big( (x^6+15 x^4-80 x^3+27 x^2-528 x+781) \sqrt{ x^4+10 x^2-96 x-71} \Big. \\ {} - \Big .(x^8 + 20 x^6 - 128 x^5 + 54 x^4 - 1408 x^3 + 3124 x^2 + 10001) \Big) + C$$
A: I have not seen anyone else mention Irresistible Integrals by Boros & Moll, so look up that book.
A: Integrals containing one or more parameters are the ideal toys for doing this sort of stuff. Besides integration and differentiation w.r.t. to the parameters, analytic continuation can yield surprising results. Thing is that when we use such methods to compute a given integral, we actually engage in such a milking exercises, albeit it with the goal of reaching the given target. 
E.g. I show here how to get to the result:
$$\int_{0}^{\infty}\frac{\log(x)}{x^2 -1}dx = \frac{\pi}{4}$$
Starting from the well known result:
$$\int_{0}^{\infty}\frac{x^{-p}}{1+x}dx = \frac{\pi}{\sin(\pi p)}$$
Differentiation w.r.t. to $p$ is an obvious step. Getting the minus 1 in the denominator when you start out with a denominator that is always nonzero requires introducing another parameter obtained by rescaling the integration variable, but this then needs to be analytically continued. 
The nice thing about this derivation is that it avoids quite some mathematical manipulations needed to do this directly via contour integration. Moving the parameter in the complex plane, moves the pole of the integrand and if we start to push the pole against the integration contour that needed to be considered for a direct derivation, then the fact that the result is analytically continued, implies that the contour would need to be deformed to make way for the approaching pole.
So, we can avoid having to even think about choosing the right contours, taking principal values etc. etc. when manipulating integrals using parameters.
