3

I have a system of ODE's and I want to check if a function of the variables is a conserved quantity of the system.

My system of ODE's is similar to this one:

x'[t_] = -b x[t] y[t]
y'[t_] =  b x[t] y[t] - d y[t]

Now, assume I want to check whether there exists constants c1 and c2 such that the following function is a conserved quantity for the system

V[t_] = c1 x[t]^2 + c2 y[t]^2

meaning, I want to find c1 and c2 such that the derivative of V w.r.t. time is zero.

I tried solving for $c_1$ and $c_2$ in the following way 

Reduce[D[V[t], t] == 0, {c1, c2}]

but this results in 

(x[t] == 0 && d == 0) || (d - b x[t] != 0 && 
c2 == (b c1 x[t])/(-d + b x[t])) || (d == 0 && b == 0 && 
x[t] != 0) || (x[t] != 0 && b == d/x[t] && d != 0 && c1 == 0) || 
y[t] == 0

which say that c1 and c2 might change with time.

How can I tell Mathematica that c1 and c2 are constants?

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
rpa
  • 287
  • 1
  • 7
  • Welcome! Thanks for taking the [tour]. It help us to help you when you write an excellent question. You can always edit if improvable. Showing due diligence, brief context, and minimal working example of code and data in formatted form is most appreciated. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. – rhermans Jul 12 '18 at 17:24
  • 1
    Thanks for the warm welcome. I look forward to contributing to the community. – rpa Jul 12 '18 at 17:33
  • 2
    For starters, don't confuse Set (=) with SetDelayed (:=) and Equal (==). Equations are defined using Equal, functions are defined with a Blank (_) pattern as an argument and Set or SetDelayed, i.e f[x_]:=1+x^2. [Edit] your question to fix that first. – rhermans Jul 12 '18 at 17:42
  • Have a read on this question and answers: What are the most common pitfalls awaiting new users? – rhermans Jul 12 '18 at 17:44
  • 1
    Is x1 supposed to be x1? Also, A trivial solution is c1 == c2== 0. – Jens Jul 12 '18 at 17:47
  • I'm referring to the x1 in the definition of V... – Jens Jul 12 '18 at 17:56
  • @Jens Sorry, yes that way a typo. The solution c1 == c2 ==0 is valid, but I am looking for a method to solve for c1 and c2 for other functions V that might be a lot more complicated. – rpa Jul 12 '18 at 18:46

2 Answers2

3

Certainly, the system of ODEs in the question has a first integral, because it is autonomous. However, the expression given in the queston for the first integral may not be correct, because differentiating it does not yield zero. Consider instead

c == -d Log[E^(-b x[t]/d) x[t]] + b y[t];

or equivalently, c == b (x[t] + y[t]) - d Log[x[t]], with c an arbitrary constant. It is a first integral, as can be seen from

Simplify[D[%, t] /. {x'[t] -> -b x[t] y[t], y'[t] -> b x[t] y[t] - d y[t]}]
(* True *)

This first integral can be derived as follows. Because t does not appear explicitly in the equations, let x be a function of y instead. Then, the two equations collapse into one.

(D[x[y], y] y'[t] == -b x[y] y) /. y'[t] -> b x[y] y - d y;
Simplify[-#/y] & /@ %
(* (d - b x[y]) x`[y] == b x[y] *)

DSolve[%, x[y], y][[1, 1]] // Simplify
(* x[y] -> -d ProductLog[-b E^((b y - C[1])/d)/d]/b *)

Equal @@ Solve[Equal @@ %, C[1]][[1, 1]] /. {C[1] -> c, x[y] -> x[t], y -> y[t]}
(* c == -d Log[E^(-b x[t]/d) x[t]] + b y[t] *)
bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • Thank you! This certainly finds the conserved quantity for the proposed SIR system. However, it relies on the insight that $dy/dx$ can be explicitly solved for using DSolve. My real system is a generalization of the system with more than two ODEs, so I'm not able to use an analogous ode like $dy/dx$ to find the analytical solution. Do you know if your procedure could be generalized to systems with more than two ODEs? – rpa Jul 13 '18 at 13:47
  • 1
    @rpa If your system of equations is autonomous (and sometimes even if it is not), then a first integral exists. However, even for just two first order autonomous equations, the resulting first integral may be a numerical function (think ParametricNDSolve with one parameter), not a symbolic function. If a numerical first integral is acceptable to you, then, yes, my approach is generalizable to more than two autonomous ODEs. There is considerable literature on this subject, largely based on Lie groups. – bbgodfrey Jul 13 '18 at 13:56
  • Thanks for the quick reply. I am interested in symbolic functions so unfortunately the numerical function would not be that useful for me. I have encountered some of the Lie group approach to first integrals along the way but according to my poor understanding of the topic, they seem to rely on knowledge about the group symmetries for the system. – rpa Jul 13 '18 at 14:12
  • 1
    Since there is only a finite number of low order Lie groups, you could try every one of them, although that could be tedious for large systems. (I used that method in my PhD thesis 50 years ago to find new solutions to Einstein's field equations.) May I suggest that you ask a new question with a three-ODE system to see what you get. Also, consider asking a question (with a different equation) on [Math.SE], where readers are likely to approach the problem from a different perspective. – bbgodfrey Jul 13 '18 at 18:27
2

The example in the question only has the trivial solution c1==c2==0 because the assumed form of V doesn't actually correspond to an integral of the motion.

To illustrate the method for finding the constants in general, let's pick a different example where at least the conserved quantity has a form similar to V with non-trivial coefficients.

I'll be using the Euler equations for the angular velocity of a freely spinning object in the body frame, where the two known conserved quantities (energy and magnitude of the angular momentum) have exactly the quadratic form given in the question. The parameters are the principal moments of inertia, i1, i2, i3. They take the place of a, b, d. The equations of motion are again similar to the ones in the question, except that there's no linear term (the original problem is actually of the Lotka-Volterra type).

After defining the derivatives and the quadratic form V, I use SolveAlways to determine the parameters. Everything that doesn't appear in the second argument of SolveAlways is considered to be a constant parameter. This is an important ingredient in finding the solution.

x'[t_] := (i2 - i3)/i1 y[t] z[t]; 
y'[t_] := (i3 - i1)/i2 x[t] z[t]; z'[t_] := (i1 - i2)/i3 x[t] y[t]

V[x_, y_, z_] := c1 x^2 + c2 y^2 + c3 z^2

sols = 
 SolveAlways[0 == D[V[x[t], y[t], z[t]], t], {x[t], y[t], z[t]}]

{{i1 -> 0, i3 -> 0}, {i1 -> 0, i3 -> i2}, {i2 -> 0, i1 -> 0}, {i2 -> 0, i3 -> 0}, {i2 -> 0, i3 -> i1}, {i2 -> i1, i3 -> 0}, {c2 -> 0, i2 -> 0}, {i2 -> i1, i3 -> i1}, {c1 -> 0, i1 -> 0}, {c2 -> c1, i2 -> i1}, {c3 -> (i3 (c2 i1^2 - c1 i2^2 - c2 i1 i3 + c1 i2 i3))/( i1 (i1 - i2) i2)}}

Inspecting the result, only the last case is of interest. The others are possible ways of getting a conserved quantities for special choices of the principal moments. The last answer is the most general.

Therefore, define V with the last choice of parameters as our conserved quantity:

Clear[c1, c2, i1, i2, i3, constOfMotion]; 
constOfMotion[c1_, c2_] = Simplify[V[x[t], y[t], z[t]] /. Last[sols]]

c1 x[t]^2 + c2 y[t]^2 + ( i3 (c2 i1 (i1 - i3) + c1 i2 (-i2 + i3)) z[t]^2)/(i1 (i1 - i2) i2)

This is the general answer for the rigid body problem. It contains two free parameters c1 and c2 because no matter what principal moments you choose there are always two constants of the motion, $I_1 x^2 + I_2 y^2 + I_3 z^2$ and $I_1^2 x^2 + I_2^2 y^2 + I_3^2 z^2$. The two parameters c1, c2 encode the weights with which these two are added to form another constant of motion.

Jens
  • 97,245
  • 7
  • 213
  • 499