Closed-form of $\int_0^1 B_n(x)\psi(x+1)\,dx$

The paper

V. Adamchik, Polygamma functions of negative order, J. Comp. Appl. Math. 100, (1998), 191-199)

contains the following statement (Proposition 3, Eq. (17)):

Let $n\in\mathbb{Z}_{\ge0}$ and $\Re z>0$, then \begin{align}\int_0^z x^n\psi\left(x\right) dx=&(-1)^{n-1}\zeta'(-n)+\frac{(-1)^n}{n+1}B_{n+1}H_n+\\ & \sum_{k=0}^n(-1)^kz^{n-k}{n \choose k}\left[\zeta'(-k,z)-\frac{B_{k+1}(z)H_k}{k+1}\right],\end{align} where $B_n$ and $B_n(z)$ are Bernoulli numbers and polynomials, $H_n$ are harmonic numbers, and $\zeta(s,z)$ denotes Hurwitz zeta function.

Specializing to $z=1$, we get $$\int_0^1 x^n\psi\left(x\right) dx=\sum_{k=0}^{n-1}(-1)^k{n \choose k}\left[\zeta'(-k)-\frac{B_{k+1}H_k}{k+1}\right].\tag{1}$$ Using (1) and the recursion relation $\psi(x+1)-\psi(x)=x^{-1}$, the integral $\int_0^1 P\left(x\right)\psi\left(x+1\right) dx$ can be computed for any polynomial $P(x)$.

Edit: There remains the task of explaining why Bernoulli polynomials provide a better basis than monomials $x^n$ (instead of linear combination of $\zeta'(-k)$ we get only one of these values). I think this could be understood after careful reading of the paper.

This is an answer for the odd case.

Proposition. Let $k=1,2,3,\ldots$. Then

$$ \int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx=(-1)^{k+1}\frac{(2k+1)!}{(2\pi)^{2k+1}}\pi \: \zeta(2k+1)-\sum_{j=0}^{2k}\!\frac{ {{2k+1}\choose j} B_j}{2k+1-j} \quad (*) $$

The proof is here, observing that $$\begin{align} \int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx & = \int_0^{1} \frac{B_{2k+1}\left(\left\{ 1/x \right\}\right)}{x}dx. \tag1 \end{align}$$ A proof of $(1)$ may be found here.

The complete answer for the question with $\,n\ge 1\,$ is the union of the answers of Start wearing purple and Olivier Oloa (the sign of the right side of the equation $(*)$ is wrong there), we get:

$$\int\limits_0^1 B_n(x)\psi(x+1)dx = (-1)^{n-1} \left(n \,\zeta’(1-n)- B_n H_{n-1}\right) + \sum\limits_{k=1}^n {\binom n k} \frac{B_{n-k}}{k}$$

A well-functioning trick is to invert $(1)$ of Start wearing purple’s answer by using: $$b_{m-1}=\sum\limits_{v=0}^{m-1} {\binom m v} a_v \enspace \Leftrightarrow \enspace a_m=\frac{1}{m+1}\sum\limits_{v=0}^m {\binom {m+1} {v+1}} B_{m-v}b_v $$

Note 1:

$$\sum\limits_{k=1}^n {\binom n k} \frac{B_{n-k}}{k}=\int\limits_0^1\frac{B_n(x)-B_n(0)}{x}dx$$

Note 2:

$$n\zeta’(1-n)- B_n H_{n-1}=\frac{d}{dt}B_t(1)|_{t=n}-H_n B_n(1)$$

$\hspace{4cm}$ Informations about $\,B_t(x)\,$ can be found here .