Consider the following equation, subject to homogeneous Neumann boundary condition. $$ u_t = \Delta u + f(u). $$ The weak formulation is as follows: $$ (u_t,w) = (\Delta u, w) + (f(u),w), \quad \forall w \in H^1(\Omega). $$ The standard finite element scheme is： $$ (u_{h,t},w_h) = (\Delta_h u_h, w_h) + (f(u_h),w), \quad \forall w_h \in V^r_h(\Omega), $$ where the $V^r_h$ denotes the space of polynomials of degree at most r. And the mass lumping method scheme is $$ (u_{h,t},w_h)_h = (\Delta_h u_h, w_h) + (f(u_h),w)_h, \forall w_h \in V^r_h(\Omega), $$

or

$$ (u_{h,t},w_h)_h = (\Delta_h u_h, w_h)_h + (f(u_h),w)_h, \forall w_h \in V^r_h(\Omega), $$

where $(\cdot, \cdot)_{h}$ signifies the use of the mass lumping method for integration. The main difference between the above two mass lumping methods is the assembly of the stiffness matrix.

Defined $R_h$：$L^2(\Omega) \to V_h(\Omega)$, and $\Pi_h: C(\bar{\Omega}) \to V^r_h$ denote the interpolation by a function in $V^r_h$ through the nodes of mass lumping integration nodes,Then we have $$ \partial_t \Pi_h u(t) - \Delta_h \Pi_h u(t) = \Pi_h f(u) + \underbrace{\Pi_h \Delta u(t) - \Delta_h \Pi_h u (t)}_{g_h(u)}. $$

$$ \partial_t \Pi_h u(t) - \Delta_h \Pi_h u(t) = \Pi_h f(u) + \underbrace{R_h \Delta u(t) - \Delta_h \Pi_h u (t)}_{g_h(u)}. $$

In the paper https://link.springer.com/article/10.1007/s10915-021-01746-y Lemma 2.4 proof that
$$ (\Pi_h \Delta u - \Delta_h \Pi_h u , w_h) \leq h^{r+1}||u||_{2r+2} ||w_h||1 $$

Consider the function $g_h(u)$ as: $$ \begin{aligned} (g_h(u),w_h) &= (R_h \Delta u - \Delta_h \Pi_h u, w_h) = (R_h \Delta u - \Delta u,w_h) + (\Delta u - \Delta_h \Pi_h u,w_h). \end{aligned} $$ Can similar results be derived? My objective is to obtain the optimal convergence rate.