Lately, I have found several interesting problems involving Harmonic numbers such as \begin{equation*}\sum_{n=1}^{\infty}\frac{H_n^{(2)}}{n^2}=\frac{7\pi^4}{360}\end{equation*} I am not familiar with computing sums involving Harmonic numbers. Is there a general strategy for tackling such problems?
How can this series be evaluated by operating on the generating function for $H_n^{(2)}$?
You may evaluate it without using generating functions.
$$S=\sum\limits_{n=1}^{\infty} \frac{H_n^{(2)}}{n^2} = \sum\limits_{n=1}^{\infty}\sum\limits_{k=1}^{n}\frac{1}{n^2k^2}$$
By changing the order of summation, you may write it as:
$$\begin{align}S &= \sum\limits_{k=1}^{\infty}\sum\limits_{n=k}^{\infty}\frac{1}{n^2k^2}\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\sum\limits_{n=k}^{\infty}\frac{1}{n^2}\right)\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\sum\limits_{n=1}^{\infty}\frac{1}{n^2} - \sum\limits_{n=1}^{k-1}\frac{1}{n^2}\right) \\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\zeta(2) - \sum\limits_{n=1}^{k-1}\frac{1}{n^2}\right)\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\zeta(2) + \frac{1}{k^2} - H_k^{(2)}\right)\\&= \zeta(2)\sum\limits_{k=1}^{\infty}\frac{1}{k^2} + \sum\limits_{k=1}^{\infty}\frac{1}{k^4} - \sum\limits_{k=1}^{\infty}\frac{H_k^{(2)}}{k^2} \\&= \zeta^2(2)+\zeta(4) - S\end{align}$$
Hence, $\displaystyle S = \frac{\zeta^2(2)+\zeta(4)}{2}$.
In general for any sequence $(a_n)_{n \ge 1}$ such that $\displaystyle \sum_{n=1}^{\infty} a_n$ converges absolutely,
We have $$\sum_{n=1}^{\infty}\left(a_n\sum_{k=1}^{n}a_k\right) = \frac{1}{2}\left(\sum_{n=1}^{\infty}a_n^2 + \left(\sum_{n=1}^{\infty}a_n\right)^2\right)$$ by using the same method as above.
Here is a general version I discovered some time ago.
Proposition : $$\sum_{r=1}^{n} \dfrac{H_{r} ^{(m)}}{r^m} = \dfrac{1}{2} \left( [H_{n}^{(m)}]^2 + H_{n}^{(2m)} \right) \quad ; \quad m \geq 1$$
Proof : Expand the summation, using the definition of $\displaystyle H_{r}^{(m)}$ , to see that $$\displaystyle \sum_{r=1}^{n} \dfrac{H_{r}^{(m)}}{r^m} = \sum_{r=1}^{n} \dfrac{1}{r^{2m}} + \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj} = H_{n}^{(2m)} + \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj} $$
And,
$\displaystyle [H_{n}^{(m)}]^2 = H_{n}^{(2m)} + 2\sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj} $
Eliminating $\displaystyle \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj}$ from the above equations, we have,
$$\displaystyle \sum_{r=1}^{n} \dfrac{H_{r}^{(m)}}{r^m} =\dfrac{1}{2} \left( [H_{n}^{(m)}]^2 + H_{n}^{(2m)} \right) \quad \square \tag{*} $$
Now put $m=2$ take limit to infinity to get the desired result.
We can also derive the result using Summation By Parts.
By using the formula in the answer to question Calculating alternating Euler sums of odd powers we have: \begin{equation} {\bf H}_2^{(2)}(t) = Li_4(t) + \frac{1}{2} [Li_2(t)]^2 +\underbrace{\int\limits_0^1 \frac{\log(\xi)}{\xi} [\log(1-t \xi)]^2 d\xi}_{I(t)} \end{equation} Now from Compute an integral containing a product of powers of logarithms. we have: \begin{equation} I(1) = \frac{1}{2}\left( -3! \zeta(4) + 2[\zeta(2)]^2\right) \end{equation} where $\Psi$ is the polygamma function. If we want to find the value at $t=-1$ it is more complicated. Here we use the identity $\log(1+\xi) = \log(1-\xi^2) - \log(1-\xi)$ and change variable in one of the integrals accordingly and use the approach developed in Calculating alternating Euler sums of odd powers. Then we have: \begin{eqnarray} I(-1) &=& \left(\frac{1}{4}-1\right) \cdot \int\limits_0^1 \frac{\log(\xi)}{\xi} \log(1-\xi)^2 d\xi - 2 \int\limits_0^1 \frac{\log(\xi)}{\xi} \log(1-\xi) \log(1+\xi) d\xi\\ I(-1)&=&\left(\frac{1}{4}-1\right) \cdot I(1) + 2 \left({\bf H}^{(1)}_3(-1) + {\bf H}_2^{(2)}(-1) +\frac{1}{2} \zeta(2)^2\right) \end{eqnarray} Bringing all together we have: \begin{eqnarray} {\bf H}^{(2)}_2(1) &=& \frac{7 \pi^4}{360}\\ {\bf H}^{(2)}_2(-1) &=& - \frac{37 \pi^4}{1440} - 2 {\bf H}^{(1)}_3(-1) \end{eqnarray}
