I want to to numerically solve the equation

\(i\partial_tU=-(\partial_x^2+\partial_y^2)U-n_2IU+n_mI^mU\)

for the field \(U(x,y,t)\), where \(I=UU^*\) and \(n_2,n_m\) are constants, m>4.

There are two integrals of motion for this equation:

\(E=\int{I}dxdy\)

\(H=\int{(\partial_\mu U \partial^\mu U-n_2I^2/2+n_mI^m/m)}dxdy\)

But the standard numerical scheme of splitting into physical factors (used for example in Spin-Glass Model Governs Laser Multiple Filamentation by Wahb Ettoumi, Jérôme Kasparian, Jean-Pierre Wolf) actually does not preserve the second integral of motion.

Where can I find in the literature a discussion of conservation properties of numerical schemes for the nonlinear Schrödinger equation? Has someone already solved the problem of preserving both conserved integrals of motion?