# Looking for an appropriate method of NDSolve for dynamical system

I am dealing with the numerical solutions of a variety of dynamical systems with integrals of motion.

As the example, let me consider the Kuramoto model, which equations of motion are

KuramotoEquations = Table[D[f[i][t], t] == g/M* Sum[Sin[f[j][t] - f[i][t]], {j,1, M}], {i, 1, M}];


where M is the number of oscillators, g is the positve constant and f[i] is the phase of oscillator. This model has a set of motion integrals (it is the rigorous fact), that is tightly related to the Mobius equation. I try to probe these integrals of motion numerically.

The Mobius transformation preserves the also-called cross-ratio. The cross ratio is defined as follows,

crossratio[f1_, f2_, f3_, f4_] := (Exp[I*f1]-Exp[I*f3])/(Exp[I*f1]-Exp[I*f4])*(Exp[I*f2] - Exp[I*f4])/(Exp[I*f2] - Exp[I*f3]);


Due to existence of the mentioned motion integrals, the whole dynamics of model can be described as the Mobius transformation that acts on the initial conditions, which I can choose randomly from the given interval,

KuramotoInitials = Table[f[i][0] == RandomReal[{0, 2*Pi}], {i, 1, M}];


In order to check that the cross-ratio is conserved, I compute it at moment t=0 for my initial conditions,

KuramotoTestCrossRatio =
crossratio[KuramotoInitialPhases[[1]], KuramotoInitialPhases[[2]],
KuramotoInitialPhases[[3]], KuramotoInitialPhases[[4]]] // Chop


Next, I perform numerical solution of the system on the interval [0, T] with pre-defined T and for a fixed constant g. I obtain the set of function f[i][t] and compute the cross-ratio as the function of time for my pre-choosen oscillators,

KuramotoCrossRatio[t_] :=
Re[crossratio[f[1][t], f[2][t], f[3][t], f[4][t]]]


From analytical treatment, I can show that KuramotoCrossRatio[t] should be simply the constant, which coincides with KuramotoCrossRatioTest. However, testing all the methods from this question, I concluded that all the methods are imperfect, despite any manipulations with WorkingPrecitions and PrecisionGoal. I have obtained the best results with implicit methods (see the fig. below).

As was noticed by @MichaelE2 , it is more reasonable to investigate the error propagation. Let me define the relative error as follows,

(KuramotoCrossRatio[t] - KuramotoTestCrossRatio)/KuramotoTestCrossRatio * 100


and denote it as $$\Delta I$$. This quantity is plotted below.

So, it seems that LinearImplicitMidpoint is quite good for this system.

As I understand, for a system with an integral of motion implicit methods and leapfrog methods are appropriate to capture the integral of motion. From my old question, I know that there is the method "Projection" by it requires to specify an explicit form of an integral of motion.

So, my question is following: consider that I have a dynamical system and I has a set of possible "candidates" for integral of motion, how should I use NDSolve in order to test my "candidates"?

If I will be completely sure that they are integrals of motion, I can simply use "Projection" method or "SymplecticPartitionedRungeKutta". However, I does not know my "candidates" are indeed integrals of motion or not.

Why numerical check is needed?

The existence of such integrals of motion was proven only for the described the simplest possible case. More complicated case corresponds to the following equation

KuramotoEquations = Table[D[f[i][t], t] == g/M* Sum[A[[i,j]] * Sin[f[j][t] - f[i][t]], {j,1, M}], {i, 1, M}];


where A[[i,j]] represents "topology" of interaction. The attempt to investigate integrals of motion is based on a conjecture that in general case there are several Mobius transformations. To verify this conjecture, one should check how many different cross-ratios exists in the system (roughly speaking).

• Can you not simply check the integral using any accurate solution produced by NDSolve? Plot[integral[t] /. sol, {t, a, b}] to see if constant? Or Plot[integral[t] - integral[a] /. sol, {t, a, b}] to see the error drift? Commented Jul 3, 2022 at 14:38
• Try FixedStep with implicit method and small stepsize
– I.M.
Commented Jul 5, 2022 at 1:50

If I understand you question correctly, you have some differential equation, and some candidates for integrals of motion, and you want to check if they indeed are integrals of motion.

Your idea is to make a numerical check. Now I have not been able to understand from your question why you insist on a numeric check. It seems to me it would be much easier to use Mathematica to check this symbolically. I will do this here for Kuramoto and cross ratio. If you do not want to make a symbolic check, and really insist on a numeric check, then please explain why.

So I take your equations but write them as replacement rules:

M = 6; (* just an example *)
KuramotoEquationsRules = {D[f[i_][t],t] :> g/M*Sum[Sin[f[j][t]-f[i][t]],{j,1,M}]};


crossratio[f1_, f2_, f3_, f4_] := (Exp[I*f1]-Exp[I*f3])/(Exp[I*f1]-Exp[I*f4])*(
Exp[I*f2]-Exp[I*f4])/(Exp[I*f2]-Exp[I*f3]);


So in this case, a candidate for an integral of motion (actually a known integral of motion) would be for example the following crossratio of the first four components:

candidate = crossratio[f[1][t],f[2][t],f[3][t],f[4][t]];


We can now symbolically check that it is conserved:

D[candidate,t] /. KuramotoEquationsRules // Simplify


which gives 0.

Perhaps you can use Mathematca to make a similar symbolic calculation for your other differential equations and your other candidates. If a symbolic calculation is not possible, and the calculation really has to be numerical for some reason, then please clarify why.

• I think it's clear the OP already knows this. I don't know why the need or desire for numerical methods needs to be clarified. It's quite clear that they are being asked for. Commented Jul 4, 2022 at 20:30
• Thanks. I agree that OP very likely knows this (at least for this equation and integral, not so sure about other equations). I am mainly trying to get OP to state more clearly what one actually wants, which is not clear to me. Currently, numerical methods appear throughout the post, but the actual question is about testing candidates for integrals of motion, not about numerical methods. If this is primarily a question about comparing various numerical schemes and their properties, then perhaps OP can formulate the question accordingly. Commented Jul 4, 2022 at 20:56
• @MichaelE2 could you please clarify what do you mean as "highgly accurate solution"? Just try yo improve WorkinPrecision and/or other built-in tools? Commented Jul 4, 2022 at 21:32
• @user293787 , very good moment! I avoided any symbolical operations. I will try to test your suggestion Commented Jul 4, 2022 at 21:34
• @ArtemAlexandrov WorkingPrecision would be the first thing I would try. Commented Jul 4, 2022 at 21:34