A combined analytical-numerical method based on a matching asymptotic algorithm is proposed for treating angular (sharp corner or wedge) singularities in the numerical solution of the Navier-Stokes equations. We adopt an asymptotic solution for the local flow around the angular points based on the Stokes flow approximation and a numerical solution for the global flow outside the singular regions using a finite-volume method. The coefficients involved in the analytical solution are iteratively updated by matching both solutions in a small region where the Stokes flow approximation holds. Moreover, an error analysis is derived for this method, which serves as a guideline for the practical implementation. The present method is applied to treat the leading-edge singularity of a semi-infinite plate. The effect of various influencing factors related to the implementation are evaluated with the help of numerical experiments. The investigation showed that the accuracy of the numerical solution for the flow around the leading edge can be significantly improved with the present method. The results of the numerical experiments support the error analysis and show the desired properties of the new algorithm, i.e. accuracy, robustness and efficiency. Based on the numerical results for the leading-edge singularity, the validity of various classical approximate models for the flow, such as the Stokes approximation, the inviscid flow model and the boundary layer theory of varying orders are examined. Although the methodology proposed was evaluated for the leading-edge problem, it is generally applicable to all kinds of angular singularities and all kinds of finite-discretization methods. © 2004 John Wiley and Sons, Ltd.