How to find Sum of Sum of Partial Harmonic Series: $\sum_{b=1}^{m} \sum_{k=1}^{b} \frac{1}{k} $? How do I find the Sum of the Sum of Finitely Many Harmonic Series: $\sum_{b=1}^{m} \sum_{k=1}^{b} \frac{1}{k} $?
According to maple it is:
$\sum_{b=1}^{m} \sum_{k=1}^{b} \frac{1}{k} =  \left(  \left( m+1 \right) ^{2}-m-4 \right)  \left( \Psi \left( m+2
 \right) +\gamma \right) + \left( - \left( m+1 \right) ^{2}+2\,m+5
 \right)  \left( \Psi \left( m+3 \right) +\gamma \right) -3
$
 A: The iterated sum actually evaluates to
$$\sum_{b = 1}^m \sum_{k = 1}^b \frac{1}k = \left(m + 1\right)H_m - m.$$
You can show by regrouping terms that for a general sum
$$\sum_{a = 1}^x \sum_{b = 1}^a f(b) = \sum_{k = 1}^x \left(x - k + 1\right) f(k).$$
So plugging in $f(b) = \frac{1}b$, we obtain
\begin{align}
\sum_{a = 1}^x \sum_{b = 1}^a f(b) &= \sum_{a = 1}^x \sum_{b = 1}^a \frac{1}b \\
&= \sum_{k = 1}^x \left(x - k + 1\right) \frac{1}k \\
&= \sum_{k = 1}^x \left(\frac{x}k - 1 + \frac{1}k\right) \\
&= \left(x + 1\right)H_x - x.
\end{align}
I'm guessing that Maple just used a weird polygamma identity somewhere.
A: Let $S_n$ denotes $\displaystyle\sum_{k=1}^n H_k$ and by applying Abel's summation:
$$\displaystyle\sum_{k=1}^n a_k b_k=A_nb_{n}+\sum_{k=1}^{n-1}A_k\left(b_k-b_{k+1}\right)\ $$ where $\ \displaystyle A_n=\sum_{i=1}^n a_i\ $ and letting $\ \displaystyle a_k=1 $ , $\ \displaystyle b_k=H_k\ $, we get,
\begin{align}
S_n&=\left(\sum_{i=1}^n1\right)H_{n}+\sum_{k=1}^{n-1}\left(\sum_{i=1}^k1\right)\left(H_k-H_{k+1}\right)\\
&=nH_{n}+\sum_{k=1}^{n-1}(k)\left(-\frac1{k+1}\right), \quad\color{blue}{\sum_{k=1}^{n-1}\frac{k}{k+1}=\sum_{k=0}^{n-1}\frac{k}{k+1}}\\
&=nH_{n}-\sum_{k=1}^n\frac{k-1}{k}\\
&=nH_{n}-\sum_{k=1}^n1+\sum_{k=1}^n\frac{1}{k}\\
&=nH_n-n+H_n
\end{align}
A: Let us use
$$H_k=- \sum_{m=1}^{k}(-1)^k\frac {{k \choose m}}{m},~~ {n \choose r}={n-1 \choose r-1}+{n-1 \choose r},~~ \sum_{k=m}^{n} {k \choose m}= {n+1 \choose m+1}.$$ Then
$$S_n=\sum_{k=1}^{n} H_k=- \sum_{k=1}^{n} \sum_{m=1}^{n} (-1)^m \frac{{n \choose m}}{m}= -\sum_{m=1}^{n} (-1)^{m} \frac{{n+1 \choose m+1}}{m}
$$ $$\Rightarrow S_n= -\sum_{m=1}^{n} (-1)^m \frac{{n \choose m}}{m}- \sum_{m=1}^{n}(-1)^m \frac{{n \choose m+1}}{m}=H_n- \sum_{m=1}^{n} (-1)^m \frac{n}{m(m+1)} {n-1 \choose m+1}$$ $$\Rightarrow S_n=H_n-n\sum_{m=0}^{n} (-1)^m \left( \frac{1}{m}-\frac{1}{m+1} \right) {n-1 \choose m}$$
$$\Rightarrow S_n= H_n-n\sum_{m=1}^{n-1} (-1)^m \frac{{n-1 \choose m}}{m}-\sum_{m=1}^{n}(-1)^m  \frac{n}{m+1} {n-1 \choose m}=H_n+nH_{n-1}-\sum_{m=1}^{n} (-1)^m {n \choose m+1}$$
$$ \Rightarrow S_n= H_n+nH_{n-1}+\sum_{p=2}^{n} (-1)^p {n \choose p}=H_n+n(H_n-\frac{1}{n})+1-n=(n+1)H_n-n.$$
A: Below I'll show three ways to utilize a generated rule for substitution, to attain the goal. One uses simplify, one uses applyrule, and one uses subs.
Firstly, Maple can convert harmonic(z) to sum form, but the reverse operation needs help. I'll construct an equation for use in substitution, from what Maple does actually compute for the value of the Sum form of harmonic(z).
restart;
convert(harmonic(b), Sum):

    Sum(b/_k1/(_k1+b),_k1 = 1 .. infinity)

value(%);

               gamma + Psi(1 + b)

So let's construct a formula to utilize,
T1 := expand(value(convert(harmonic(z),Sum)) = harmonic(z));

                    1                      
     T1 := Psi(z) + - + gamma = harmonic(z)
                    z                      

Performing only the inner summation,
Q1 := Sum(sum(1/k, k=1..b), b=1..m);

       Sum(gamma+Psi(1+b),b = 1 .. m)

We can now simplify with T1 as a so-called "side-relation", after substituting z for b in Q1,
simplify(subs(b=z,  combine(expand(Q1))), {T1});

       Sum(harmonic(z),z = 1 .. m)

At which point Maple does the following,
value(%);

      harmonic(m + 1) (m + 1) - m - 1

Alternatively we can put T1 in a slightly different form, isolating for Psi(z),
T2 := isolate(T1, Psi(z));

                                          1
     T2 := Psi(z) = harmonic(z) - gamma - -
                                          z

Now take the double summation,
temp := sum(sum(1/k, k=1..b), b=1..m):

Q2 := simplify(expand(temp),size);

         / 2    \                           2                
         \m  + m/ Psi(m) + 1 + (gamma - 1) m  + (gamma + 1) m
   Q2 := ----------------------------------------------------
                                  m                          

And we can construct a rule from T2, for use with the applyrule command. This does not entail substituting m for z.
R := subs(Psi(z)=Psi(z::anything), T2):

normal(applyrule(R, Q2));

        harmonic(m) m + harmonic(m) - m

There are other ways to use T2 to substitute for Psi(z). By first substituting m for z in T2 we can use subs directly.
normal(subs(subs(z=m, T2), Q2));

        harmonic(m) m + harmonic(m) - m

