A novel method is proposed to calculate the area of vulnerability (AOV) to voltage sags, which can reflect the relationship between the position of fault occurrence and voltage sags at sensitive buses. Previous methods based on interpolation approach are mainly based on the assumption that the fault voltage magnitude curve has a unimodal shape, which does not hold in all cases. Moreover, the accuracy of polynomial interpolation is largely affected by the selection strategy of sample points, and inappropriate initial points will lead to divergence. To overcome these drawbacks, an effective approach which does not rely on sample results is proposed, where the fault voltage magnitude curve could be any arbitrary shapes. The AOV is derived by solving the parametric fault voltage equations, where Galerkin method is used to identify the polynomial approximation coefficients. With the order of polynomial approximation increased, an iteration scheme is applied, where the results in each iteration are selected as the initial points of the next iteration. Case studies show that the proposed method can provide very accurate results with high efficiency.