A generalized Simpson's rule In solving this problem, I developed a generalized Simpson's rule
$$\int_a^bf(t)dg(t)=\begin{bmatrix}f(a),\;f(\frac{a+b}{2}),\;f(b)\end{bmatrix}\left(\begin{array}{rrr}
\displaystyle-\frac{1}{2} & \displaystyle\frac{2}{3} & \displaystyle-\frac{1}{6}\\\displaystyle-\frac{2}{3} & 0\, & \displaystyle\frac{2}{3}\\
\displaystyle\frac{1}{6} & \displaystyle-\frac{2}{3} & \displaystyle\frac{1}{2}
\end{array}\right)\begin{bmatrix}g(a)\\g(\frac{a+b}{2})\\g(b)\end{bmatrix}+\mathcal{O}(|a-b|^5).$$
assuming $f(t)$ and $g(t)$ are sufficiently differentiable. This formula is exact for $f(t)=t^m$ and $g(t)=t^n$ with $m+n\leq 4$, where $m,n$ are nonnegative integers. The formula is invariant under integration by parts, i.e., it exactly reproduces
$$\int_a^bf(t)dg(t)+\int_a^bg(t)df(t)=f(b)g(b)-f(a)g(a).$$
In case $g(t)=t$ is linear, the method reduces to the standard Simpson's rule. Is this formula seen in the literature? Please show me a reference.
 A: What you have here is a special case of a big class of differential operators called 'Summation-by-Parts' (SBP) operators. They are defined as follows:
A matrix $D$ is an SBP operator of order $s$ if $D$ can be decomposed as $D = P^{-1}Q$ such that the following properties hold:
1) $D\mathbf{x}^m = m\mathbf{x}^{m-1}, \quad m=0, 1, \dots, q$,
2) $P = P^\top > 0$,
3) $Q + Q^\top = B = \text{diag}(-1, 0, \dots, 0, 1)$.
For two vectors $\mathbf{f}$ and $\mathbf{g}$, we may define the discrete inner product $(\mathbf{f},\mathbf{g})_P = \mathbf{f}^\top P \mathbf{g}$. Then, since $D \mathbf{g} \approx \mathbf{g}_x$, the identity
$$
(\mathbf{f},D \mathbf{g})_P = \mathbf{f}^\top Q \mathbf{g} = \mathbf{f}^\top (B-Q^\top) \mathbf{g} = f_n g_n-f_0g_0 - \mathbf{f}^\top Q^\top \mathbf{g} = f_n g_n-f_0g_0 - (D\mathbf{f},\mathbf{g})_P
$$
is the discrete version of the integration by parts rule
$$
(f,g_x) = \int_a^b fg_x dx = f(b)g(b)-f(a)g(a) - \int_a^b f_xgdx = f(b)g(b)-f(a)g(a) - (f_x,g),
$$
as you noted.
While SBP methods can be constructed for almost any type of differential operators (finite difference, finite volume, finite element etc.), the particular one you have derived originates within the context of spectral collocation methods. Let $P_n(x)$ be the Legendre polynomial of order $n$. Further, let $(x_0, \dots, x_n)$ be the roots of $P'_n(x)(1-x^2)$ (i.e. the Legendre-Gauss-Lobatto nodes) and let $L_j(x), j=0, \dots, n$ be the Lagrange interpolation polynomials associated with these points. Then $P$ is the diagonal matrix whose non-zero elements are given by
$$
P_{j,j} = \int_{-1}^1 L_j(x) dx
$$
and the elements of $D$ are given by
$$
D_{i,j} = L'_j(x_i).
$$
The matrix you have derived (with $n=2$) is $Q$ and is simply given by $PD$, or explicitly
$$
Q_{i,j} = \int_{-1}^1 L_i(x) L'_j(x) dx.
$$
Of course, the method may be generalised to any $n \geq 1$.
For a detailed study of these methods, see Carpenter and Gottlieb (1995), Spectral methods on arbitrary grids.
