at the four points $t^{n-j}$, $0\leq j\leq 3$. (If $0\leq n\leq 2$, we use the same linear combination, putting $U^j=U^0$ if $j<0$.)

The integrals involving the local terms are computed in general using the 5-point Gauss-Legendre quadrature rule in each spatial interval. The inner product $(P_N\mathcal{G}{u_h}_x,\chi_x)$ involving the nonlocal term is computed as the inner product $(I_N\mathcal{G}{u_h}_x,\chi_x)$ where the Fourier interpolant $I_N$ is defined as $$\begin{aligned}
\label{E312}
I_N v(x)=\sum_{k=-N/2}^{N/2-1}\hat{v}_ke^{{\rm i}kx},\end{aligned}$$ where by $\hat{v}_k$ we denote the discrete Fourier coefficients of $v$, computed by the Fast Fourier Transform. The inner product $(\cdot,\cdot)$ is approximated by the trapezoidal quadrature rule, which is very accurate for periodic functions.

In the sequel, we shall use the fully discrete scheme described above with the $C^{2}$ cubic splines ($r=4$) as the finite element subspace $S_{h}$. We shall refer to this method as the [*hybrid scheme/method*]{}.

We checked numerically the orders of convergence of the hybrid scheme as follows. Due to lack of analytical formulas for solutions of the Benjamin equation we considered the nonhomogeneous equation $$\label{E313}
u_t+uu_x+\mathcal{G} u_{xx}+\frac{1}{2}u_{xxx}=f(x,t),\quad (x,t)\in [-1,1]\times [0,T],$$ with periodic boundary conditions and $$f(x,t)=e^t\left(\sin(\pi x)+\frac{\pi}{2}e^t\sin(2\pi x)+\left(\pi^2-\frac{\pi^3}{2}\right)\cos(\pi x) \right).$$ The specific equation has a solution $u(x,t)=e^t \sin(\pi x)$. We solved it numerically up to 