AbstractHigher‐order p and hp finite element algorithms for modeling the flow of New‐tonian fluids have been developed. Numerical computations in three test geometries were carried out to compare the performance of these methods to that of conventional h‐type finite elements. Computations demonstrated that for solutions which are analytic everywhere in the domain of interest, p extensions exhibit an exponential rate of convergence as opposed to the characteristic linear convergence of h‐type finite element procedures. Moreover, for solutions which are analytic everywhere in the domain except at a set of points, hp extensions with geometric mesh refinement based on effective error indicators (individually, the error in the energy norm, local and global deviations from mass and momentum conservation) also give rise to an exponential convergence rate. This rapid convergence of p and hp finite element formulations enables the achievement of the solution of linear Stokes problems with an order of magnitude reduction in degrees of freedom. In nonlinear problems, these computational savings are appreciated by a factor approximately equal to the number of iterations required to obtain converged solutions.