Numerical study of fluid interfaces is a difficult task due to the presence of high frequency numerical instabilities. Small perturbations even at the round-off error level may experience rapid growth. This makes it very difficult to distinguish the numerical instability from the physical one. Here, we perform a careful numerical stability analysis for both the spatial and time discretization. We found that there is a compatibility condition between the numerical discretizations of the singular integral operators and of the Lagrangian derivative operator. Violation of this compatibility condition will lead to numerical instability. We completely eliminate the numerical instability by enforcing this discrete compatibility condition. The resulting scheme is shown to be stable and convergent in both two and three dimensions. The improved method enables us to perform a careful numerical study of the stabilizing effect of surface tension for fluid interfaces. Several interesting phenomena have been observed. Numerical results will be presented.