How find this sum $\sum\limits_{n=0}^{\infty}\frac{1}{(3n+1)(3n+2)(3n+3)}$
Your try is excellent... Note that $f(0)=f'(0)=f''(0)=0$ hence $$ f(1)=\int_0^1f'(x)\,\mathrm dx=\int_0^1\int_0^xf''(y)\,\mathrm dy\,\mathrm dx=\int_0^1\int_0^x\int_0^yf'''(z)\,\mathrm dz\,\mathrm dy\,\mathrm dx, $$ that is, $$ 2f(1)=\int_0^1(1-z)^2f'''(z)\,\mathrm dz=\int_0^1\frac{1-z}{1+z+z^2}\,\mathrm dz. $$ The rest is routine. The change of variable $2z+1=\sqrt3t$ yields $$ 2f(1)=\int_{1/\sqrt3}^\sqrt3\frac{\sqrt3-t}{1+t^2}\,\mathrm dt=\left[\sqrt3\arctan t-\frac12\log(1+t^2)\right]_{1/\sqrt3}^\sqrt3. $$ Note that $\arctan\sqrt3=\pi/3$ and $\arctan1/\sqrt3=\pi/6$, hence $$ 2f(1)=\sqrt3\cdot\left(\frac\pi3-\frac\pi6\right)-\frac12\log4+\frac12\log\frac43, $$ that is, $$ f(1)=\frac14\left[\frac\pi{\sqrt3}-\log3\right]\approx0.1788. $$ Second method: The rational fraction is such that $$ \frac2{(3n+1)(3n+2)(3n+3)}=\frac1{3n+1}-\frac2{3n+2}+\frac1{3n+3}, $$ hence $$ 2f(x)=x^2g_1(x)-2xg_2(x)+g_3(x),\qquad g_k(x)=\sum_{n\geqslant0}\frac{x^{3n+k}}{3n+k}. $$ Thus, for each $k$, $$ g'_k(x)=\sum_{n\geqslant0}x^{3n+k-1}=\frac{x^{k-1}}{1-x^3}. $$ Since $g'_k(0)=0$ for every $k\geqslant1$, this yields $$ 2f(x)=x^2\int_0^x\frac{1}{1-t^3}\mathrm dt-2x\int_0^x\frac{t}{1-t^3}\mathrm dt+\int_0^x\frac{t^2}{1-t^3}\mathrm dt. $$ The change of variable $t=xu$ yields $$ 2f(x)=x^3\int_0^1\frac{1}{1-x^3u^3}(1-2u+u^2)\mathrm du, $$ that is, $$ 2f(x)=x^3\int_0^1\frac{1-u}{1-xu}\frac{1-u}{1+xu+x^2u^2}\mathrm du. $$ When $x\to1$, one obtains once again $$ 2f(1)=\int_0^x\frac{1-u}{1+u+u^2}\mathrm du. $$ More generally, for every integer $k\geqslant2$, $$ \sum_{n\geqslant0}\frac{(k-1)!}{(kn+1)(kn+2)\cdots(kn+k)}=\int_0^1\frac{(1-u)^{k-2}}{1+u+\cdots+u^{k-1}}\,\mathrm du. $$
We all know that sums of the form $\displaystyle{\sum\frac1{(n+a)(n+b)}}$ or $\displaystyle{\sum\frac1{(n+a)(n+b)(n+c)}}$ etc., are telescopic in nature, and solving them is trivial: for integer values of a, b, and c, that is ! But what if a, b, and c are not integers? What then? Suddenly, things aren't so simple and trivial anymore, and the formerly banal and deceitfully tame-looking problem takes on whole new twist, gaining entirely unexpected dimensions of depth, meaning, and insight. The key, when dealing with such surprising turns of events, is to just take a step back, and try to parse or rephrase the older, well-worn solution in terms which might prove useful or relevant to circumventing the obstacles raised by the newly-encountered situation.
For instance, it is by no means difficult to show that $\displaystyle{\sum_{n=1}^\infty\frac1{(n+a)(n+b)}}=\frac{H_a-H_b}{a-b}$ , or that $$\sum_{n=1}^\infty\frac1{(n+a)(n+b)(n+c)}=\frac{(a-b)\cdot H_c\ +\ (b-c)\cdot H_a\ +\ (c-a)\cdot H_b}{(a-b)(b-c)(c-a)}$$ for natural values of a, b, c. Indeed, the user Random Variable has already proved it on this thread, though I now realize that my initial formula is slightly mistaken, in the sense that either the sum should begin at $0$, with the $\gamma$ present, or at $1$, with the $\gamma$ absent. Now, the whole question becomes how to redefine $H_m$ , so as to be able to extend its meaning to non-natural arguments. For natural arguments, we have the well-known formula $H_m=\displaystyle\sum_{k=1}^m\frac1k$ . Now, let us take the simple function $f_k(x)=\displaystyle\frac{x^k}k$ , and notice that $f_k'(x)=x^{k-1}$ . And since we know that $\displaystyle{\sum_{k=1}^mx^{k-1}=\frac{1-x^m}{1-x}}$ , this finally allows us to conclude that $\displaystyle{H_m=\int_0^1\frac{1-x^m}{1-x}dx}$ , which, unlike the previous formula, can easily be extended to non-natural arguments as well. This has already been done by Euler two and a half centuries ago, so it's hardly new territory. Now, in our case, $\{a,b,c\}=\left\{\frac13,\frac23,1\right\}$ , and the value of $H_1$ is $1$, so all that's left to do is to compute $H_\frac13$ and $H_\frac23$ using the above formula, since $$\sum_{n=0}^\infty\frac1{(3n+1)(3n+2)(3n+3)}=\frac16+\frac1{27}\cdot\sum_{n=1}^\infty\frac1{\left(n+\frac13\right)\left(n+\frac23\right)(n+1)}$$
$$H_\frac13=\int_0^1\frac{1-\sqrt[3]x}{1-x}dx=\int_0^1\frac{1-t}{1-t^3}d\left(t^3\right)=\int_0^1\frac1{1+t+t^2}d\left(t^3\right)=\int_{\sqrt[3]0}^{\sqrt[3]1}\frac{3\,t^2}{1+t+t^2}dt=$$
$$=3\int_0^1\left(1-\frac{1+t}{1+t+t^2}\right)dt=3\,\bigg[\int_0^11\cdot dt-\tfrac12\int_0^1\frac{2+2t}{1+t+t^2}dt\bigg]=$$
$$=3\,\bigg[t|_0^1-\tfrac12\bigg(\int_0^1\frac1{1+t+t^2}dt+\int_0^1\frac{1+2t}{1+t+t^2}dt\bigg)\bigg]=$$
$$=3\,\bigg[1-\tfrac12\int_0^1\frac1{\left(t+\frac12\right)^2+\frac34}dt-\tfrac12\cdot\ln\left(1+t+t^2\right)_0^1\bigg]=$$
$$=3\,\bigg[1-\tfrac12\cdot\frac1{\sqrt3/2}\cdot\arctan\bigg(\frac{t+\frac12}{\sqrt3\big/2}\bigg)_0^1-\frac{\ln3}2\bigg]=3-\frac\pi{2\sqrt3}-\tfrac32\ln3.$$
Similarly for $\displaystyle{H_{\frac23}=\tfrac32+\frac\pi{2\sqrt3}-\tfrac32\ln3}$ . Then, by substituting these values back into the original formula above, we finally arrive at the desired result $\displaystyle{I=\frac{\pi\sqrt3-3\ln3}{12}}$ .
$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle #1 \right\rangle} \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace} \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack} \newcommand{\dd}{{\rm d}} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\ic}{{\rm i}} \newcommand{\imp}{\Longrightarrow} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\pars}[1]{\left( #1 \right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert #1 \right\vert} \newcommand{\yy}{\Longleftrightarrow}$ $\ds{ I_{N} \equiv \sum_{n = 0}^{N}{1 \over \pars{3n + 1}\pars{3n + 2}\pars{3n + 3}}\,, \qquad I = I_{\infty} =\ {\large ?}}$
\begin{align} I_{N}&\equiv {1 \over 6}\sum_{n = 0}^{N}{1 \over n + 1/3} - {1 \over 3}\sum_{n = 0}^{N}{1 \over n + 2/3} + {1 \over 6}\sum_{n = 0}^{N}{1 \over n + 1} \\[3mm]&= {1 \over 6}\sum_{n = 0}^{N}\pars{{1 \over n + 1/3} - {1 \over n + 2/3}} + {1 \over 6}\sum_{n = 0}^{N}\pars{{1 \over n + 1} - {1 \over n + 2/3}} \\[3mm]&= {1 \over 18}\sum_{n = 0}^{N}{1 \over \pars{n + 1/3}\pars{n + 2/3}} - {1 \over 18}\sum_{n = 0}^{N}{1 \over \pars{n + 1}\pars{n + 2/3}} \end{align} We'll use the well known identity: $\ds{\sum_{n = 0}^{\infty}{1 \over \pars{n + z_{0}}\pars{n + z_{1}}} ={\Psi\pars{z_{0}} - \Psi\pars{z_{1}} \over z_{0} - z_{1}}}$ where $\Psi\pars{z}$ is the $\it\mbox{digamma function}$
\begin{align} I & = I_{\infty} = {1 \over 18}\,{\Psi\pars{1/3} - \Psi\pars{2/3} \over 1/3 - 2/3} - {1 \over 18}\,{\Psi\pars{1} - \Psi\pars{2/3} \over 1 - 2/3} \\[3mm]&= {1 \over 6}\,\bracks{\Psi\pars{2 \over 3} - \Psi\pars{1 \over 3}} - {1 \over 6}\,\bracks{\Psi\pars{1} - \Psi\pars{2 \over 3}} = -\,{1 \over 6}\,\Psi\pars{1} - {1 \over 6}\,\Psi\pars{1 \over 3} + {1 \over 3}\,\Psi\pars{2 \over 3} \end{align} However, \begin{align} \Psi\pars{1} &= -\gamma\quad\mbox{where}\ \gamma\ {\it\mbox{is the Euler-Mascheroni constant}}\ \pars{~\gamma = 0.5772\ldots~} \\[3mm] \Psi\pars{1 \over 3}&=-\gamma - {\root{3} \over 6}\,\pi - {3 \over 2}\,\ln\pars{3} \\[3mm] \Psi\pars{2 \over 3}&=-\gamma + {\root{3} \over 6}\,\pi - {3 \over 2}\,\ln\pars{3} \\ I & = I_{\infty} = -\,{1 \over 6}\pars{-\gamma} - {1 \over 6}\bracks{-\gamma - {\root{3} \over 6}\,\pi - {3 \over 2}\,\ln\pars{3}} + {1 \over 3}\bracks{-\gamma + {\root{3} \over 6}\,\pi - {3 \over 2}\,\ln\pars{3}} \\[3mm]&= {\root{3} \over 12}\,\pi - {1 \over 4}\,\ln\pars{3} \end{align}
$$\color{#0000ff}{\large \sum_{n = 0}^{N}{1 \over \pars{3n + 1}\pars{3n + 2}\pars{3n + 3}} = {\root{3} \over 12}\,\pi - {1 \over 4}\,\ln\pars{3}} \approx 0.1788 $$