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), $$$$ (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), $$
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), $$$$ (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), $$
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 $P_h$:$L^2(\Omega) \to V_h(\Omega)$ denote the $L_2$-projection onto the space of $V_h^r$ , 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. What I want to do is construct the error equation $e_h = u_h - \Pi_h u$,Then we have then $$ \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{P_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 Consider the function $g_h(u)$ asin the second equation: $$ (g_h(u),w_h) = (P_h \Delta u - \Delta_h \Pi_h u, w_h) = (P_h \Delta u - \Delta u,w_h) + (\Delta u - \Delta_h \Pi_h u,w_h). $$
$$ =(\Delta_h R_h u - \Delta u, w_h) + (\Delta u - \Delta_h \Pi_h u,w_h) $$
$$ =(\Delta u - \Delta_h \Pi_h u,w_h). $$ where the $R_h$ is the Ritz projection. These two schemes seem to be equivalent, but when I use $(\cdot,\cdot)_h$ to integrate the stiffness matrix, the convergence order will be reduced in some cases. Did I go wrong somewhere?