Numerical approximation of the Ericksen-Leslie system with variable density is considered in this paper. The spherical constraint condition of the orientation field is preserved by using polar coordinates to reformulate the system. The equivalent new system is computationally cheaper because the vector function of the orientation field is replaced by a scalar function. An iteration penalty method is applied to construct a numerical scheme so that stability is improved. We first prove that the scheme is unique solvable and unconditionally stable in energy. Then we show that this scheme is of first-order convergence rate by rigorous error estimation. Finally, some numerical simulations are performed to illustrate the accuracy and effectiveness of the scheme.