In this paper, we develop a gradient recovery based linear (GRBL) finite element method (FEM) and a Hessian recovery based linear FEM for second‐order elliptic equations in non‐divergence form. The elliptic equation is casted into a symmetric non‐divergence weak formulation, in which second‐order derivatives of the unknown function are involved. We use gradient and Hessian recovery operators to calculate the second‐order derivatives of linear finite element approximations. Although, thanks to low degrees of freedom of linear elements, the implementation of the proposed schemes is easy and straightforward, the performances of the methods are competitive. The unique solvability and the H2$$ {H}^2 $$ seminorm error estimate of the GRBL scheme are rigorously proved. Optimal error estimates in both the L2$$ {L}^2 $$ norm and the H1$$ {H}^1 $$ seminorm have been proved when the coefficient is diagonal, which have been confirmed by numerical experiments. Superconvergence in errors has also been observed. Moreover, our methods can handle computational domains with curved boundaries without loss of accuracy from approximation of boundaries. Finally, the proposed numerical methods have been successfully applied to solve fully nonlinear Monge–Ampère equations.