In this paper, a comparison of weak-form Galerkin and least-squares finite element models of Timoshenko beam theory with the von Kármán strains is presented. Computational characteristics of the two models and the influence of the polynomial orders used on the relative accuracies of the two models are discussed. The degree of approximation functions used varied from linear to the 5th order. In the linear analysis, numerical results of beam bending under different types of boundary conditions are presented along with exact solutions to investigate the degree of shear locking in the newly developed mixed finite element models. In the nonlinear analysis, convergences of nonlinear finite element solutions of newly developed mixed finite element models are presented along with those of existing traditional model to compare the performance.