FUN with f̶̶l̶̶a̶̶g̶̶s̶ Newton Cotes Quadrature formula and Bernoulli polynomials of the second kind I was told to phrase my question in a more exciting way when I asked it last time. The following is a preliminary consideration. If you don't need it, just scroll down to START HERE. Here we go then:
Imagine you have a really cool integral,
\begin{align}
    I = \int_a^b f(x) \ dx ,
\end{align}
which you wanted to solve numerically. And what do you take for it: A quadrature formula. And one of the most popular quadrature formulas are the Newton-Cotes formulas. Don't be put off, this just means that the supports for my quadrature formula are at a fixed distance from each other. It applies to an interval $[a,b]$:
\begin{align*}
            x_k = a + k \cdot h, \hspace{15pt} h = \frac{b-a}{n}, \hspace{15pt} k = 0,1,\ldots,n
\end{align*}
The quadrature error $E_n[f] := I_n - I$ of a Newton-Cotes formula $n$-th degree to $n$ support points,
\begin{align*}
    I = \int_a^b f(x) \, d x, \hspace{30pt} I_n = I_n\left(f,[a,b]\right) = \sum_{k=0}^n \, w_kf(x_k)\,,   
\end{align*}
for a $(n+1)$ times continuously differentiable function $f(x)$ has the representation
\begin{align*}
    E_n[f] &= \frac{f^{(n+1)}(\xi)}{(n+1)!} \int_a^b \omega(x) \,d x, \hspace{30pt} \omega(x) = \prod_{k=0}^n\,(x-x_k) \\
    &\leq \frac{1}{(n+1)!} \,\max_{\xi\in[a,b]}\left|f^{(n+1)}(\xi)\right|\,\int_a^b \omega(x) \,d x
\end{align*}
BUT ALL OF THIS DOESN'T GIVE A SHIT NOW. Because I want to specify the error exactly. It follows:
\begin{align*}
    \int_a^b \omega(x) \,d x &= \int_a^b (x-x_0)(x-x_1)\cdots(x-x_n)\,d x \\
    &= \int_a^b (x-x_0)(x-x_0 - h)(x- x_0 - 2h)\cdots(x-x_0-nh)\,d x \\
    &= \int_a^b h^{n+1}\,\left(\frac{x-x_0}{h}\right)\left(\frac{x-x_0}{h} - 1\right)\cdots\left(\frac{x-x_0}{h}-n\right)\,d x, \hspace{30pt} \color{red}{u = \frac{x-x_0}{h}} \\
    &= h^{n+2} \int_{\frac{a-x_0}{h}}^{\frac{b-x_0}{h}} u(u-1)\cdots(u-n)\,d u \\
    &= h^{n+2}\,(n+1)!\,\psi_{n+2}(u)\,\Big|_{u_1=\frac{a-x_0}{h}}^{u_2=\frac{b-x_0}{h}} \\
    &= h^{n+2}\,(n+1)!\,\left[\psi_{n+2}(u_2)-\psi_{n+2}(u_1)\right]
\end{align*}
And please do not run away in panic now. These $\psi_n(u)$ here are Bernoulli polynomials of the second kind, certain weird polynomials like Legendre polynomials, which have great properties that I want to use, see Wikipedia.
\begin{align}
   (n+1)! \cdot \psi_{n+2}(u) = \int u(u-1)\cdots(u-n)\ d u
\end{align}
For the limits $u_1 = 0$, $u_2 = n$ we find the following expression:
\begin{align*}
    \int_a^b \omega(x) \, d x &= h^{n+2}\,(n+1)!\,\left[\psi_{n+2}(n) - \psi_{n+2}(0)\right] \\
    &= h^{n+2}\,(n+1)!\,\left[-|G_{n+2}| - G_{n+2}\right] \color{white}{\frac{1}{2}} \\
    &= h^{n+2}\,(n+1)!\,\left[(-1)^n-1\right]\,|G_{n+2}| \color{white}{\frac{1}{2}}
\end{align*}
where $G_n$ are the Gregory coefficients defined as follows:
\begin{align*}
    \frac{x}{\ln(1+x)} = 1 + \frac{1}{2}x - \frac{1}{12}x^2 + \frac{1}{24}x^3 - \ldots = 1 + \sum_{n=1}^\infty\,G_n\,x^n
\end{align*}
This means that for even $n$ we cannot find an expression. Thus, for even $n$, we must use a different approximation.  For this we double the term $x_{n/2} = \frac{a+b}{2}$. So we calculate in such a way that we consider an additional grid point and thus a double zero for the interpolation polynomial.
\begin{align*}
    \int_a^b \omega(x)\,d x &= h^{n+3} \int_{0}^{n} u(u-1)\cdots(u-\tfrac{n}{2})^2\cdots(u-n)\,d u \\
    &= h^{n+3} \bigg[\underbrace{\int_0^{n} u\cdot u(u-1)\cdots(u-n)\,d u}_{\color{red}{=\,I}} - \frac{n}{2}\underbrace{\int_0^{n} u(u-1)\cdots(u-n)\,d u}_{\color{red}{=\,0}}\bigg]
\end{align*}

START HERE:
And now we have finally reached the point where I need your help. Your superhuman knowledge of mathematics is now required. Unfortunately, there was probably still no mathematician who has set up a suitable relation for me, which I can use to solve this integral:
\begin{align}
I = \int_0^{n} u\cdot u(u-1)\cdots(u-n)\,d u
\end{align}
I want to use this:
\begin{align}
   (n+1)! \cdot \psi_{n+2}(u) = \int u(u-1)\cdots(u-n)\ d u
\end{align}
I had the following ideas:

*

*Integration by parts:

\begin{align}
 I &= u\cdot (n+1)!\,\psi_{n+2}(u)\,\Big|_0^{n} - (n+1)! \int_0^{n} \psi_{n+2}(u)\,d u \\
    &= (n+1)!\left[n\cdot\psi_{n+2}(n) - \int_0^n \psi_{n+2}(u)\,d u\right] \\
    &= -(n+1)!\left[n\cdot|G_{n+2}| + \int_0^n \psi_{n+2}(u)\,d u\right]
\end{align}
But what is now
\begin{align}
\int_0^n \psi_{n+2}(u)\,d u \ ?
\end{align}


*Consider the case $n+1$
\begin{align}
 (n+2)!\cdot\psi_{n+3}(u) &= \int u(x-1)(u-2)\cdots(u-n)(u-(n+1))\ du \\
 &=\underbrace{\int u \cdot u(u-1)(u-2)\cdots(u-n) \ du}_{\color{red}{=\,I(u)}} - (n+1) \underbrace{\int u(u-1)(u-2)\cdots(u-n) \ du}_{\color{red}{=\,(n+1)!\,\cdot\,\psi_{n+2}(u)}} \\
&= I(u) - (n+1)\cdot(n+1)!\,\cdot\,\psi_{n+2}(u)
\end{align}
Thus applies:
\begin{align}
 I(u) &= (n+1)\cdot(n+1)!\,\cdot\,\psi_{n+2}(u) + (n+2)! \cdot \psi_{n+3}(u) \color{white}{\frac{1}{2}} \\ \\
I &= I(u) \ |_0^n = I(n) - I(0) \\
&= (n+1)! \cdot (n+1) \,\underbrace{\left[\psi_{n+2}(n) - \psi_{n+2}(0)\right]}_{\color{red}{=\,0; \text{ $n$ is even}}} \color{white}{\frac{1}{2}} \\
&\hspace{15pt} + (n+2)! \left[\psi_{n+3}(n) - \psi_{n+3}(0)\right] \color{white}{\frac{1}{2}} \\
&= (n+2)! \left[\psi_{n+3}(n) - \psi_{n+3}(0)\right]
\end{align}
But what is now
\begin{align}
\psi_{n+3}(n) \ ?
\end{align}
 A: According to Wikipedia:
\begin{align}
    \psi_{2k}(k-1+y) = \psi_{2k}(k-1-y) 
\end{align}
If we now choose $y = k - 2$, it follows:
\begin{align}
    &\psi_{2k}(k-1+(k-2)) = \psi_{2k}(k-1-(k-2)) \color{white}{\frac{1}{2}}\\
    & \Longrightarrow \hspace{10pt} \psi_{2k}(2k-3) = \psi_{2k}(1) \color{white}{\frac{1}{2}} \\
\end{align}
If we consider an even $n$ as described in my problem, i.e. $n = 2k$, the following holds:
\begin{align}
    &\psi_{n}(n-3) = \psi_{n}(1) =  G_{n-1} + G_n  \\
&\psi_{n+3}(n) = \psi_{n+3}(1) =  G_{n+2} + G_{n+3} 
\end{align}
We get for the integral:
\begin{align}
I &= (n+2)! \cdot \left[\psi_{n+3}(n) - \psi_{n+3}(0)\right] \color{white}{\frac{1}{2}}\\
&= (n+2)! \cdot \left[G_{n+2} + G_{n+3} - G_{n+3}\right] \color{white}{\frac{1}{2}} \\
&= (n+2)! \cdot G_{n+2}
\end{align}
Thus, for Newton-Cotes formulas:
\begin{align*}
    E_n[f] = \left\{\begin{array}{ll} \displaystyle -f^{(n+1)}(\xi)\,h^{n+2}\cdot 2|G_{n+2}|, & n \text{ ungerade} \\ \displaystyle -f^{(n+2)}(\xi)\,h^{n+3}\cdot G_{n+2}, & n \text{ ungerade} 
    \end{array}\right.
\end{align*}
