Given a set of ODE's $ \dot{z} = f(z) $ (or discrete time $ z_{t+1} = f(z_t) $), is there a way to numerically find constants of motion?

For $ f(z_t) \approx M z_t $, diagonalizing the matrix $ M $ and finding eigenvalues equal to 1 works, but this only finds linear combinations of $ z $ that are constant. I would like to find things like radius, energy, and momentum.

I am specifically interested in constants of motion as I integrate Lie groups.


Just to clarify, a constant of motion, at least for the sake of this questions, is defined as,

\begin{align} f: \mathbb{R}^n \rightarrow \mathbb{R} \\ \text{s.t.}~~\frac{d}{dt} f(z(t)) = 0 \end{align}

So, for a given trajectory $ z(t) $, we would expect perhaps multiple constants of motion $ f_i(z(t)) $, but I would like to find constants of motion $ f_i $ such that they are constant for all possible trajectories.

\begin{align} \frac{d}{dt} f_i(z(t; z_0^{(k)})) = 0 ~~ \forall z_0^{(k)} \in \mathbb{R}^n \end{align}

  • $\begingroup$ Do you know beforehand how many constants of motion are you computing? $\endgroup$
    – nicoguaro
    Commented Aug 23, 2017 at 17:51
  • $\begingroup$ @nicoguaro, No. If I am integrating groups, then it should be $ \text{dim}(\mathbb{R}^n/G) $ $\endgroup$
    – user23933
    Commented Aug 23, 2017 at 18:02
  • $\begingroup$ So, are you hoping to find a method that tell you how many independent invariants do you have and also tell you the invariants? $\endgroup$
    – nicoguaro
    Commented Aug 24, 2017 at 19:38
  • $\begingroup$ Ideally yes. Practically, I would be happy with one that gets me the "first" N invariants and then I just guess. $\endgroup$
    – user23933
    Commented Aug 24, 2017 at 19:48

3 Answers 3


This seems like a difficult problem to solve in general, but here's one idea.

Suppose you look for all constants of motion of the form $$ \sum_{i=1}^m \lambda_i \varphi_i(t), $$ where $\lambda_i$ are some unknown coefficients, and $\{\varphi_i\}_i$ is some set of basis functions that can be computed by numerical integration of the ODE. For example, $\varphi_i$'s could be all monomials up to a certain degree, composed of $z,\dot z,\ddot z,\ldots$.

Then the condition that $\sum_i \lambda_i \varphi_i(t)=\mathrm{const}$ could be approximated as the matrix equation $$ \Phi \lambda = 0, \qquad \text{where }\Phi_{ij} = \varphi_j(t_i) - \varphi_j(0), $$ where $\lambda$ is a vector unknown, $\Phi$ is an $n\times m$ matrix, and $t_1,\ldots,t_n$ are some arbitrary convenient points in time ($n>m$).

This way you could find all possible coefficients $\lambda$ that generate constants of motion by asking for the numerical kernel of $\Phi$. If the constants of motions are polynomials in $z$ and its derivatives, this would work, but the matrix could be very large due to how many different monomials you can have. This does require you to guess in advance what the constants of motion might look like, I don't know if that's realistic for you.


If you can write down the problem using generic constants then this becomes a parameter estimation problem. The DifferentialEquations.jl docs explains a bunch of methods in more detail, but I'll focus on a simple and standard method exemplified here. That example is about finding biological constants (birth and death rates in Lotka-Volterra) but the idea is the same as finding constants of motion: convert the problem into an ODE with parameters and optimize the parameter values.

The inner loop is essentially:

  1. Solve the problem with parameters $p$.
  2. Compute a cost function.
  3. Have the optimization routine update the parameters $p$ to hopefully better parameters.
  4. Go back to 1 if your cost is too high.

If you have some data for what the path of the differential equation should be, then you can define a cost function that is the l2 difference between the data and the numerical solution at each point. Minimizing this cost is the same as maximizing the Normal log-likelihood, and results in finding the parameters of the ODE which best fit the data.

But this technique doesn't require data. You can define any loss function on the solution which makes sense, and have an optimizer fit it. For example, you might only know that you want the maximum value of $u$ at times $t = 1,2,3,4$ to be $1$. So you can make the cost function be

$$C(u) = \sum_{t=1,2,3,4} \Vert u(t) - 1 \Vert $$

and have the optimizer search for parameters which best satisfy this constraint. Of course, you get better results the better you define your cost function to match what you really want, though more complicated cost functions will put more of a strain on the optimizer.

Basic optimization methods here include Levenberg-Marquardt and you'll see some tutorials using that, but you really shouldn't be using that because it's very sensitive to finding local minima and these parameter estimation problems are very tough for local optimizers. Instead, you should be using a good global optimization library like NLopt to get decent results (it'll still be difficult!)

  • $\begingroup$ Once you find one constant of motion, how do you find the others? Also, is there a way to make the constants scale invariant? For example, a 3D vector undergoing 2D rotations should have $ r, \phi $ as constants. First, I'd like to find both, and furthermore find a mapping (not just the values) so that I can map large amounts of vectors to their constants. $\endgroup$
    – user23933
    Commented Aug 23, 2017 at 15:26
  • $\begingroup$ For the scale thing, I think I could just replace 1 with a free parameter per data sample and heavily regularize the other parameters to norm 1. Finding the other constants is still hard $\endgroup$
    – user23933
    Commented Aug 23, 2017 at 15:30
  • $\begingroup$ You should have the optimizer optmize the vector of all parameters together. Well, there are some cases where you may want to single out some to do independently first, or iteratively optimize some, but in general you optimize them together because the diffeq depends on all of them together. $\endgroup$ Commented Aug 23, 2017 at 15:48
  • $\begingroup$ Can you illustrate how this would work for $\dot z=Mz$? On cursory reading, it seems to me like you're answering a slightly different question, but I could be wrong. $\endgroup$
    – Kirill
    Commented Aug 23, 2017 at 16:06
  • $\begingroup$ After re-reading, I agree with Kirill. I think you are trying to find constants in the differential equation. I am trying to find constants of motion. A constant of motion is a function of the state that does not vary with time. $ \frac{d}{dt} f(z(t)) = 0 $ but depends strictly on the initial conditions $ f(z_i(t)) = f(z_i(0)) = \lambda_i $. I am not looking for the parameters of the ODE. $\endgroup$
    – user23933
    Commented Aug 23, 2017 at 16:20

Still hoping for a better answer, but here is something I have thought of (inspired by other what Christ Rackauckas said), but don't particularly like.

Define a parameterized "constant of motion" $ f_i(z(t); \theta) $ and then optimize on $ \theta $ like,

\begin{align} \mathcal{L}(\theta_i, \vec{\lambda}) &= \frac{1}{2}\big(\|\theta\|^2 - 1\big) + \sum_{t, k} \frac{1}{2} \big( f_i(z^{(k)}(t); \theta_i) - \lambda_{ki} \big)^2 \end{align}

We leave $ \lambda_{ki} $ "free" so that each constant for each trajectory, can still scale, but is subject to be a function of the state at each time step.

Training for the next constant of motion could be done by adding a loss term to account for the last constant of motion. \begin{align} \sum_{t, k} \sum_{j=1}^{i-1} \big( f_i(z^{(k)}(t); \theta_i) - f_j(z^{(k)}(t); \theta_j) \big)^2 \end{align}

I don't like this for the obvious reason that it's a horrible optimization problem to solve.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.