Nuclear symmetry energy $E_{rm{sym}}(rho)$ at density $rho$ is normally expanded or simply parameterized as a function of $chi=(rho-rho_0)/3rho_0$ in the form of $E_{rm{sym}}(rho)approx S+Lchi+2^{-1}K_{rm{sym}}chi^2+6^{-1}J_{rm{sym}}chi^3+cdots$ using its magnitude $S$, slope $L $, curvature $K_{rm{sym}}$ and skewness $J_{rm{sym}}$ at the saturation density $rho_0$ of nuclear matter. Much progress has been made in recent years in constraining especially the $S$ and $L$ parameters using various terrestrial experiments and astrophysical observations. However, such kind of expansions/parameterizations do not converge at supra-saturation densities where $chi$ is not small enough, hindering an accurate determination of high-density $E_{rm{sym}}(rho)$ even if its characteristic parameters at $rho_0$ are all well determined by experiments/observations. By expanding the $E_{rm{sym}}(rho)$ in terms of a properly chosen auxiliary function $Pi_{rm{sym}}(chi,Theta_{rm{sym}})$ with a parameter $Theta_{rm{sym}}$ fixed accurately by an experimental $E_{rm{sym}}(rho_{rm{r}})$ value at a reference density $rho_{rm{r}}$, we show that the shortcomings of the $chi$-expansion can be completely removed or significantly reduced in determining the high-density behavior of $E_{rm{sym}}(rho)$. In particular, using two significantly different auxiliary functions, we show that the new approach effectively incorporates higher $chi$-order contributions and converges to the same $E_{rm{sym}}(rho)$ much faster than the conventional $chi$-expansion at densities $lesssim3rho_0$. Several quantitative demonstrations using Monte Carlo simulations are given.