Let $u_h$ be the finite element solution of a fourth order equation (like biharmonic equation), using polynomial degree two.If the convergence rate of $u_h$ is $2$, what is the convergence rate of the second order derivatives of $u_h$? Is it possible they do not converge?
Convergence of the second derivative of the finite element solution
finite element;error estimation
Yes, for several reasons.First, it is instructive to look at how one usually proves convergence rates. In the standard setting, where your discrete approximation $u_h\in V_h$ is in the same space $V$ as the true solution $u$, the error $\|u-u_h\|_V$ in the natural space $V$ (where you have continuity and coercivity of the bilinear form) is bounded by the interpolation error of $u$ by functions in $V_h$ (this is Ca's lemma). So all you have to do is to get an estimate of the interpolation error $\|u- I_h u\|_V$. If $V=H^l$, the space of $l$ times weakly differentiable functions, and $V_h$ is the space of piecewise polynomials of degree $k$ on a given mesh, and your exact solution is in $H^{k+1}$ for $k+1\geq l$, then you can show that$$ \| u- I_h u \|_{H^l} \leq C h^{k+1-l} |u|_{H^{k+1}},\tag{1}$$where $|u|_{H^{s}}$ is the $H^s$ seminorm of $u$, i.e., the $L^2$ norm of all partial derivatives of order $s$, and $C>0$ is a constant independent of $h$ and depending only on the geometrical properties of the mesh. (The proof works by first considering a reference element, where the error can be bounded by a constant, and then transforming the norms to each element and summing; the power of $h$ in the estimate comes purely from how the different $H^s$ norms scale under affine transformations.) This means that the smoothness of the true solution limits the polynomial degree you can use to get better convergence. (You can always use higher degrees, of course, but you won't get better accuracy from them.) You can use different methods (e.g., Aubin-Nitsche lemma) to trade a lower order error norm for higher powers of $h$, i.e., an estimate of the form $(1)$ for the error $\|u-u_h\|_{H^l}$ usually holds also for smaller values of $l$. Now to your specific case. The natural regularity of the biharmonic equation is $H^2$, since the weak form involves second derivatives. You are also using piecewise quadratic polynomials. This means $l=k=2$, such that you would need to have the true solution in $H^3$, which is more than the natural regularity, and hence does not necessarily hold for the test problem you are solving. (But $u\in H^2$ is enough to get from $(1)$ a convergence order of $\mathcal{O}(h^2)$ in $H^0=L^2$, which is what you observe. You should also get $\mathcal{O}(h)$ in $H^1$, i.e., for the first derivatives, which is a useful test.)Rule of thumb: You lose one order of convergence for each derivative, and you gain one order for each polynomial degree, up to a point determined by the regularity of the exact solution you are trying to approximate.However, you are not in the standard setting: Functions in $H^2$ are  continuously differentiable across element boundaries (apply the usual argument to the first derivatives), which is not the case for the standard piecewise polynomial spaces (which are only continuous across boundaries). Hence, $(1)$ doesn't apply in this case. Intuitively, since you are not enforcing the correct continuity across element boundaries, you might get the function values right (enough), but not the derivatives -- and especially not the second derivatives. (Here's where the test with the first derivatives comes in handy.) Instead, you also need an interior penalty method to enforce the required continuity of $H^2$ functions across element boundaries (such as discontinuous Galerkin methods; an alternative would be mixed methods). But how to do this would be a new question.