In the dynamical systems approach to describing turbulent or otherwise chaotic flows, an important quantity is the Lyapunov exponents and vectors that characterize the strange attractor of the flow. In particular, knowledge of the Lyapunov exponents and vectors will help identify perturbations that the system is most sensitive to, and quantify the dimension of the attractor. However, reliably computing these quantities requires robust numerical algorithms. While several perturbation-based techniques are available in literature, their application to commonly used turbulent flow solvers as well the numerical convergence properties have not been studied in detail. The goal of this work is to examine the convergence properties of the Lyapunov exponents when computed using a low Mach number solver, for a series of progressively complex flow problems. In particular, a manufactured solutions approach is devised using the Orr-Sommerfeld (OS) perturbation theory to extract Lyapunov exponents and vectors from OS eigen solutions. Overall the spatial convergence rates of Lyapunov exponents follow the truncation order for the discretization scheme. However, temporal convergence rates were found to be only a weak function of timestep used. Additionally, for some configurations, the convergence properties are found to depend on the Lyapunov exponent itself. These results indicate that convergence properties do not follow universal rates, and require careful analysis for specific configurations considered.