Integration of a function approximated by a nth order polynomial I've been playing with Simpson's rule and a thing came up to my mind. The rectangular rule is a 0th order polynomial approximation of integration. The trapezoidal rule is 1st. Simpson's rule is 2nd. Then what about nth? I've been working on this issue more than 2 weeks and couldn't get the answer. How could I solve this?
What I did is below.
approach 1
A nth order polynomial function passing $n+1$ points $(x_i, y_i) (0 \leq i\leq n)$ can be written by Lagrange polynomial as 
$$ g(x) = \sum_{j=0}^{n} y_i \displaystyle \prod_{i=0 \atop i \neq j}^{n} \frac{x-x_i}{x_j - x_i}$$
So, if $y=f(x)$, $g(x)$ is a polynomial approximation of $f(x)$. So, the integration of $f(x)$ can be written as follow.
$$ \int_{x_0}^{x_n} f(x) dx \simeq \int_{x_0}^{x_n} g(x) dx = \int_{x_0}^{x_n} \sum_{j=0}^{n} f(x_i) \displaystyle \prod_{i=0 \atop i \neq j}^{n} \frac{x-x_i}{x_j - x_i} dx \\
= \sum_{j=0}^{n} f(x_i) \displaystyle \prod_{i=0 \atop i \neq j}^{n} (x_j - x_i)^{-1} \int_{x_0}^{x_n} \prod_{i=0 \atop i \neq j}^{n} (x-x_i) dx
$$
Here, we suppose $x_i = x_0 + hi$ where $h$ is constant, $z = \frac{x - x_0}{h}$ and $z_i = \frac{x_i - x_0}{h} = i$ and the equation can be simplified.
$$\sum_{j=0}^{n} f(x_i) \displaystyle \prod_{i=0 \atop i \neq j}^{n} (x_j - x_i)^{-1} \int_{x_0}^{x_n} \prod_{i=0 \atop i \neq j}^{n} (x-x_i) dx \\
=  \sum_{j=0}^{n} f(x_i) \displaystyle \prod_{i=0 \atop i \neq j}^{n} (j - i)^{-1} h \int_{0}^{n} \prod_{i=0 \atop i \neq j}^{n} (z-i) dz$$
$f(x_i)$ was left intact intentionally. Solving equation boils down to solving $\int_{0}^{n} \displaystyle \prod_{i=0 \atop i \neq j}^{n} (z-i) dz$ which is quite difficult. 
$S(z, j) := \displaystyle \prod_{i=0 \atop i \neq j}^{n} (z-i)$ can be written as follows.
$$S(z, j) = \displaystyle \prod_{i=0 \atop i \neq j}^{n} (z-i)
= \sum_{i=0}^{n} a_i z^i$$
$a_i$ is constant. Notice the upper value of $\sum$ is $n$ for the number of $z$s in $\prod$ is $(n - 0 + 1) - 1 = n$ because $i \neq j$. 
It is trivial $S(z, j) = 0$ when $z = 0, 1, ... , n$ but $\neq j$. So the following is true.
$$\sum_{i=0}^{n} a_i z^i = 0 \,(z = 0, 1, ... , n \cap z \neq j)$$
The integration can be calculated as below.
$$\int_{0}^{n} \displaystyle \prod_{i=0 \atop i \neq j}^{n} (z-i) dz = 
\int_{0}^{n} S(z, j) dz = \int_{0}^{n} \sum_{i=0}^{n} a_i z^i dz \\
= \big[\sum_{i=0}^{n} a_i \frac{z^{i+1}}{i+1}\big]_{0}^{n} \\
= \sum_{i=0}^{n} a_i \frac{n^{i+1}}{i+1}$$
To solve this, evaluating $a_i$ is required. $a_i$ can be evaluated as below.
$$\sum_{i=0}^{n} a_i z^i = 0 \,(z = 0, 1, ... , n \cap z \neq j) \\
\Leftrightarrow
\displaystyle 
  {\begin{pmatrix}
    1 & 0 & 0 & \cdots & 0 \\
    1 & 1 & 1 & \cdots & 1 \\
    1 & 2 & 2^{2} & \cdots & 2^{n} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    1 & n & n^{2} & \cdots & n^{n} 
  \end{pmatrix}}
{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\\vdots \\a_{n}\end{pmatrix}}
= {\begin{pmatrix}b_{0}\\b_{1}\\b_{2}\\\vdots \\b_{n}\end{pmatrix}} \\
(b_i =   
\begin{cases}
  S(j,j)\,\,(i = j) \\
  0 \,\,(i \neq j)
\end{cases}
)
$$
As $a_n = 1$ and $a_0 = \prod_{i=0 \atop i \neq j}^{n} -i$ the equation can be written as follows.
$$\sum_{i=0}^{n} a_i z^i = 0 \,(z = 0, 1, ... , n \cap z \neq j) \\
\Leftrightarrow
\sum_{i=0}^{n-1} a_i z^i = -z^n \,(z = 1, 2, ... , n \cap z \neq j) \\
\Leftrightarrow
\displaystyle 
  {\begin{pmatrix}
    1 & 1 & 1 & \cdots & 1 \\
    1 & 2 & 2^{2} & \cdots & 2^{n-1} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    1 & n & n^{2} & \cdots & n^{n-1} 
  \end{pmatrix}}
{\begin{pmatrix}a_{1}\\a_{2}\\\vdots \\a_{n-1}\end{pmatrix}}
= {\begin{pmatrix}b_{1}\\b_{2}\\\vdots \\b_{n-1}\end{pmatrix}} \\
(b_i =   
\begin{cases}
  S(j,j) - z^j\,\,(i = j) \\
  -z^i \,\,(i \neq j)
\end{cases}
)
$$
The leftmost matrix is Vandermonde's matrix and multiplying its inverse matrix to the left of both sides gives $a_i$. However, the inverse matrix is quite complicated and I find it pretty hard to calculate (Inverse of Vandermonde's matrix).
approach 2
I guessed $i \neq j$ makes the problem difficult. So I went for deleting it.
$$ \int_{0}^{n} \displaystyle \prod_{i=0 \atop i \neq j}^{n} (z-i) dz =
\int_{0}^{n} \frac{\displaystyle \prod_{i=0}^{n} z-i}{z-j} dz$$
The limit of $\frac{\prod_{i=0}^{n} z-i}{z-j}$ when $z$ approaches $j$ is given by L'hospital's rule, which is finite.
$$ \frac{(\prod_{i=0}^{n} z-i)'}{(z-j)'}|_{z=j} = (\prod_{i=0}^{n} z-i)'|_{z=j}$$
But after hours of thinking, I couldn't come up with how to solve $\int_{0}^{n} \frac{\displaystyle \prod_{i=0}^{n} z-i}{z-j} dz$. Two $z$s at both the enumerator and denominator. Integration by parts didn't work.
I've read the following articles and still didn't make it to find the answer.


*

*Riemann sum

*Newton–Cotes formulas

*Lagrange polynomial

*Vandermonde matrix

*Gaussian quadrature
How do I solve this?
 A: That is a very interesting question. If you're still interested here are some things you can look into.
The sequence for the n-th order is given in the OEIS, namely at A002177, A002178, A002179, .... You are looking for the Cotesian numbers.
The whole triangle (this is probably the one you want, with the values calculated up to 100-th degree polynomials) is given at A100642. I believe that that is what you are looking for.
In the above sequences there is a link to a scanned copy of W.M. Johnson's paper "On Cotesian numbers: their history, computation and values to n=20". This will provide you with some explanation on how to calculate them.

Update (1 july 2022): I conducted the full derivation of the closed form expression for the Cotesian numbers in my bachelor thesis "Methods for reducing error in approximations of the Rayleigh integral" in appendix B.
If we let $s(n, k)$ denote the signed Stirling numbers of the first kind, the closed form expression of the Cotesian numbers is proven to be
$$c_i = \frac{1}{(n-1)!} \binom{n}{i} \sum_{j=0}^n \sum_{m=0}^{n-j} i^m n^j \frac{(-1)^{i+n} s(n+1, j+m+1)}{j+1}\,.$$

Update (4 october 2022): The specific pages of the derivation can now also be found on my github page.
