If $S$ is a N-a-K spline with knots $x_1, \dotsc, x_4$ then it satisfies the spline conditions: twelve equations in twelve unknowns. (Twelve coefficients, six equations to prescribe values at the knots and six more to force continuity of derivatives up to third order at $x_2$ and $x_3$.) Since $p_1, p_2, p_3$ fit up to third order in all (two) inner knots, it follows that $p_1 = p_2 = p_3$. So $S$ and $P$ are both Lagrange interpolation polynomials and therefore $S=P$.
Instead of $4$ pairs of points, let's look at the case with $3$ pairs of points first.
- Let $H = \big\{ (x,y,h) \in \mathbb{R}^3 : h \ge 0 \big\}$ be the upper half-space.
- For any $p = (x,y,h) \in H$, let $\tilde{p} = (x,y,0)$ be its projection on the $xy$-plane.
Given any $3$ points $p_1, p_2, p_3\in H$ whose projections $\tilde{p}_1$, $\tilde{p}_2$, $\tilde{p}_3$ bounding a triangle in the $xy$-plane in counter-clockwise orientation,
consider the convex hull formed by following $6$ vertices $p_1, p_2, p_3, \tilde{p}_1, \tilde{p}_2, \tilde{p}_3$. The resulting polyhedron is a truncated triangular prism.
To compute its volume, we split it into $3$ tetrahedra:
$$\tilde{p}_1, \tilde{p}_2, \tilde{p}_3, p_3;\quad
\tilde{p}_1, \tilde{p}_2, p_3, p_2;\quad\text{ and }\quad
\tilde{p}_1, p_2, p_3, p_1$$
We find
$$\begin{align}
\verb/Volume/ &=
\frac16
\begin{vmatrix}
1 & x_1 & y_1 & 0\\
1 & x_2 & y_2 & 0\\
1 & x_3 & y_3 & 0\\
1 & x_3 & y_3 & h_3\\
\end{vmatrix}
+
\frac16
\begin{vmatrix}
1 & x_1 & y_1 & 0\\
1 & x_2 & y_2 & 0\\
1 & x_3 & y_3 & h_3\\
1 & x_2 & y_2 & h_2\\
\end{vmatrix}
+
\frac16
\begin{vmatrix}
1 & x_1 & y_1 & 0\\
1 & x_2 & y_2 & h_2\\
1 & x_3 & y_3 & h_3\\
1 & x_1 & y_1 & h_1\\
\end{vmatrix}\\
\\
&= \frac{h_1+h_2+h_3}{6}\begin{vmatrix}
1 & x_1 & y_1\\
1 & x_2 & y_2\\
1 & x_3 & y_3\\
\end{vmatrix}
\end{align}
$$
For the truncated quadrilateral prism at hand, we can split it into two
truncated triangular prisms.
Up to a sign, the volume you seek will be equal to:
$$\verb/Volume/ =
\frac{h_1+h_2+h_3}{6}\begin{vmatrix}1 & x_1 & y_1\\1 & x_2 & y_2\\1 & x_3 & y_3\end{vmatrix}
+
\frac{h_1+h_3+h_4}{6}\begin{vmatrix}1 & x_1 & y_1\\1 & x_3 & y_3\\1 & x_4 & y_4\end{vmatrix}
$$
Best Answer
While Acheille Hui has outlined a good strategy I would like to continue with your approach.
First of all, I can confirm that your polynomial passes through the four given values, hence you have the correct expression.
You can always change the interval of integration from, say, $[a,b]$ to $[-1,1]$ by a linear transformation, specifically \begin{equation} x = \phi(t) = a \frac{1-x}{2} + b \frac{1+x}{2}, \end{equation} where I have chosen the peculiar representation to stress the fact that \begin{equation} \phi(-1) = a, \quad \phi(1) = b \end{equation} We also see that \begin{equation} \phi'(t) = \frac{b-a}{2}. \end{equation} It follows that for any continuous function $f : [a,b] \rightarrow \mathbb{R}$ we can write \begin{equation} \int_a^b f(x)dx = \int_{-1}^1 f(\phi(t)) \phi'(t) dt \end{equation}
In your case we could write \begin{equation} x_1 = a, \quad x_2 = a + h, \quad x_3 = a + 2h, \quad x_4 = a + 3h = b \end{equation} so that \begin{equation} \phi(t) = a \frac{1-x}{2} + (a+3h) \frac{1+x}{2} = a + 3h \frac{1+x}{2} \end{equation} and \begin{equation} \phi(-1) = x_1, \quad \phi\left(-\frac{1}{3}\right) = x_2, \quad \phi\left(\frac{1}{3}\right) = x_2, \quad \phi\left(1\right) = x_3 \end{equation} The substitution will allow you to move your integral to the interval $[-1,1]$. At this point, I would recommend reducing the polynomials to standard form by carrying all the multiplications.
You mention Simpson's rule. It is difficult to apply it in this context, because the number of natural sub-intervals is an odd number (specifically 3), rather than an even number. The simplest way that I know to rediscover Simpson's rule is to seek numbers $a$, $b$, and $c$ such the quadrature rule \begin{equation} \int_{-h}^h f(x) dx = a f(-h) + b f(0) + c f(h) \end{equation} is exact for all polynomials of degree less than or equal to $2$. By using \begin{equation} f(x) = 1, \quad f(x) = x, \quad \text{and}\quad f(x) = 2 \end{equation} you will assemble a linear system of dimension $3$ for the unknowns $a$, $b$ and $c$. Once you have solved it you will be able to verify that your new quadrature rule also integrates $f(x) = x^3$ correctly.