The research work reported in this dissertation is aimed to develop efficient and stable numerical schemes in order to obtain accurate numerical solution for viscoelastic fluid flows within the spectral element context. The present research consists in the transformation of a large class of differential constitutive models into an equation where the main variable is the logarithm of the conformation tensor or a quantity related to it in a simple way. Particular cases cover the Oldroyd-B fluid and the FENE-P model. Applying matrix logarithm formulation in the framework of the spectral element method is a new type of approach that according to our knowledge no one has implemented before. The reformulation of the classical constitutive equation using a new variable namely the logarithmic formulation, enforces the eigenvalues of the conformation tensor to remain positive for all steps of the simulation. However, satisfying the symmetric positive definiteness of the conformation tensor during the simulation is the necessary condition for stability; but definitely, it is not the sufficient condition to reach meaningful results. The main effort of this research is devoted to introduce a new algorithm in order to overcome the drawback of direct reformulating the classical constitutive equation to the logarithmic one. To evaluate the capability of the extended matrix logarithm formulation, comprehensive studies have been done based on the linear stability analysis to show the influence of this method on the resulting eigenvalue spectra and explain its success to tackle high Weissenberg numbers. With this new method one can treat high Weissenberg number flows at values of practical interest. One of the worst obstacles for numerical simulation of viscoelastic fluids is the presence of spurious modes during the simulation. At high Weissenberg number, many schemes suffer from instabilities and numerical convergence may not be attainable. This is often attributed to the presence of solution singularities due to the geometry, the dominant non-linear terms in the constitutive equations, or the change of type of the underlying mixed-form differential system. Refining the mesh proved to be not very helpful. In this study, to understand more deeply the mechanism of instability generation a comprehensive study about the growth of spurious modes with time evolution, mesh refinement, boundary conditions and Weissenberg number or any other affected parameters has been performed. Then to get rid of these spurious modes the filter based stabilization of spectral element methods proposed by Boyd was applied with success.