In general, solutions to the Laplacian equation enjoy relatively high smoothness. However, they can exhibit singular behaviors at domain corners or points where boundary conditions change type. In this talk, I will focus on the mixed Dirichlet-Neumann boundary conditions for Laplacian equation, and discuss how singularities in this case adversely affect the accuracy and convergence rates of standard numerical methods. Then, starting from the celebrated Weierstrass theorem about polynomial approximation, I will describe the approximation theory related to the so-called Muntz polynomials, which can be viewed as a generalization of usual polynomials. Additionally, I will illustrate the idea of Muntz-Galerkin methods, and show that how they can overcome the difficulties to achieving high order accuracy for the problems with singularities.