SUMMARY The existence of pores, cracks and cleavage in rocks results in significant non-linear elastic phenomena. One important non-linear elastic characteristic is the deviation of the stress–strain curve from the linear path predicted by Hooke's law. To provide a more accurate description of the non-linear elastic characteristics of rocks and to characterize the propagation of non-linear elastic waves, we introduce the Preisach–Mayergoyz space model. This model effectively captures the non-linear mesoscopic elasticity of rocks, allowing us to observe the stress–strain and modulus–stress relationships under different stress protocols. Additionally, we analyse the discrete memory characteristics of rocks subjected to cyclic loading. Based on the Preisach–Mayergoyz space model, we develop a new non-linear elastic constitutive relationship in the form of an exponential function. The new constitutive relationship is validated through copropagating acousto-elastic testing, and the experimental result is highly consistent with the data predicted by the theoretical non-linear elastic constitutive relationship. By combining the new non-linear elastic constitutive relationship with the strain–displacement formula and the differential equation of motion, we derive the non-linear elastic wave equation. We numerically solve the non-linear elastic wave equation with the finite difference method and observe two important deformations during the propagation of non-linear elastic waves: amplitude attenuation and dispersion. We also observe wave front discontinuities and uneven energy distribution in the 2D wavefield snapshot, which are different from those of linear elastic waves. We qualitatively explain these special manifestations of non-linear elastic wave propagation.