  Nuclear symmetry energy $E_{\rm{sym}}(\rho)$ at density $\rho$ is normally
expanded or simply parameterized as a function of $\chi=(\rho-\rho_0)/3\rho_0$
in the form of $E_{\rm{sym}}(\rho)\approx
S+L\chi+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 $\lesssim3\rho_0$. Several quantitative
demonstrations using Monte Carlo simulations are given.
