The spatial discontinuity of physical parameters at an abrupt interface may increase numerical errors when solving partial differential equations. Rather than generating boundary-adapted meshes for objects with complicated geometry in the finite-element method, the subpixel smoothing (SPS) replaces discontinuous parameters inside square elements that are bisected by interfaces in, for example, the finite-difference (FD) method, with homogeneous counterparts and matches physical boundary conditions therein. In this work, we apply the idea of SPS to the eight-band effective-mass Luttinger-Kohn (LK) and Burt-Foreman (BF) Hamiltonians of semiconductor nanostructures. Two smoothing approaches are proposed. One stems from eliminations of the first-order perturbation in energy, and the other is an application of the Hellmann-Feynman (HF) theorem. We employ the FD method to numerically solve the eigenvalue problem corresponding to the multiband Schrodinger’s equation for circular quantum wires (QWRs). The eigen-energies and envelope (wave) functions for valence and conduction states in III-V circular QWRs are examined. We find that while the procedure of perturbation theory seems to have the better accuracy than that of HF theorem, the errors of both schemes are considerably lower than that without smoothing or with direct but unjustified averages of parameters. On the other hand, even in the presence of SPS, the numerical results for the LK Hamiltonian of nanostructures could still contain nonphysical spurious solutions with extremely localized states near heterostructure interfaces. The proper operator ordering embedded in the BF Hamiltonian mitigates this problem. The proposed approaches may improve numerical accuracies and reduce computational cost for the modeling of nanostructures in optoelectronic devices.