Based on the Gross–Pitaevskii equation model that incorporates higher-order nonlinear interaction, we studied sonic black hole and sonic horizon formation involved in nonlinear dynamical evolution for Bose–Einstein condensates with higher-order nonlinear interaction. On the basis of the modified variational method and the scenario where the system starts dynamic evolution from ground state, we derived the typical system distribution width function, which is analytically formulated as periodic oscillation solution and monotonically damped variational oscillation solution under different parametric settings. We also calculated the criteria formula for sonic black hole horizon formation with regard to the two evolution modes: oscillation mode and monotonically decay mode, pictorially demonstrating the time interval of sonic horizon appearance. The theoretical results obtained here can be used to guide relevant experimental studies of sonic horizon and sonic black hole formation for Bose–Einstein condensates incorporating the higher-order nonlinear interaction effects.