I am trying to reproduce the results of a paper A new test for chaos in deterministic systems by Georg A. Gottwald and Ian Melbourne. This paper talks of a simple 1-0 method to determine whether a dynamical system is chaotic, or not. It's short and simple, so I'll summarise the method below.

The code I wrote is below at the end. My $K$ value is always around 6 or 7 but it should be either near 0 or 1. No matter what I try I do not see where I made a mistake. I did as the paper said. Please help!

Assume you have the following dynamical system

enter image description here

What you need to do is calculate the integrals below, where $\phi(\textbf{x}(s))$ is any observable of the system. The procedure is independent of the choice of $\phi(\textbf{x}(s))$ and the constant $c>0$. The paper chose $\phi(\textbf{x}(t))=x_1(t)+x_2(t)$ and $c=1.7$

enter image description here

And, finally using $p(t)$ calculate

enter image description here

According to the paper, $K$ in that limit must be either $0$ (for nonchaotic systems) or $1$ (for chaotic systems). And, that concludes the test.

Now, the dynamical system that the paper uses is the forced van de Pol oscillator, whose equations of motion could be written as below. Paper chose $a=d=5$, and let $\omega$ vary from $2.462$ to $2.466$ in steps of size $0.00001$. Also, by converting the above two equations for $\theta$ and $p$ to differential equations, we have 4 differential equations which can be numerically solved.

enter image description here

x10 = 1; x20 = 1; time = 1000;
 c = 1.7; d = 5; a = 5; \[Omega] = 2.43;
inis = {x1[0] == 1, x2[0] == 1, \[Theta][0] == 0, p[0] == 0};
eqns = {
   x1'[t] == x2[t],
   x2'[t] == -d (x1[t]^2 - 1) x2[t] - x1[t] + a Cos[\[Omega] t],
   \[Theta]'[t] == c + (x1[t] + x2[t]) - (x10 + x20),
   p'[t] == (x1[t] + x2[t]) Cos[\[Theta][t]] - (x10 + x20)
{X1, X2, Th, P} = 
  NDSolveValue[Join[eqns, inis], {x1, x2, \[Theta], p}, {t, 0, time}];
K[t_, T_] := 
 Log[ 1/T NIntegrate[(P[t + \[Tau]] - P[\[Tau]])^2, {\[Tau], 0, T}]]/
K[Exp[10], Exp[10]]
  • $\begingroup$ Criterium K they proposed is too tricky, since it needs an integration time of $T=2\times 10^6$. It can be compare to $T=10^3$ you have used. I can recommend you to remove x10, x20 from the code and replace K with k since K is the symbol occupied by the system. Then add Method -> {"StiffnessSwitching", "NonstiffTest" -> False} to NDSolve and try time= 2 10^6. $\endgroup$ Nov 17, 2020 at 11:23
  • $\begingroup$ @AlexTrounev I tried adding the method of integration but it's not much of a difference. When I set time=2 10^6 it takes so long to process that it is practically impossible to plot the type of graph that the paper has plotted of K vs ω. Did you try to run the code with time= 2 10^6, was it fast? $\endgroup$ Nov 17, 2020 at 12:55
  • $\begingroup$ @AlexTrounev it finally produced an output after running for about 7 minutes or a little more but it's practically then same as when I had time=1000, that is K is around 6. I don't understand why it would give 6 when the maths I used is exactly what the paper used. $\endgroup$ Nov 17, 2020 at 13:42
  • $\begingroup$ No it is not fast. Also this method of chaotic and nonchaotic solutions determination looks very strait and not connected with numerical method in use. $\endgroup$ Nov 17, 2020 at 13:46
  • $\begingroup$ Now use second restriction underlined in the paper that $t<<T$ and use Log[M[t]+1]/Log[t] as they recommended. $\endgroup$ Nov 17, 2020 at 13:51

1 Answer 1


We can reproduce one test from the paper using next code

x10 = -1.654211168152969`; x20 = 0.6666053281765396`; time = 500000;
c = 1.7; d = 5; a = 5; \[Omega] = 2.463;
inis = {x1[0] == x10, x2[0] == x20, \[Theta][0] == 0, p[0] == 0};
eqns = {x1'[t] == x2[t], 
   x2'[t] == -d (x1[t]^2 - 1) x2[t] - x1[t] + 
     a Cos[\[Omega] t], \[Theta]'[t] == c + (x1[t] + x2[t]), 
   p'[t] == (x1[t] + x2[t]) Cos[\[Theta][t]]};
{X1, X2, Th, P} = 
  NDSolveValue[Join[eqns, inis], {x1, x2, \[Theta], p}, {t, 0, time}, 
   Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];

Kk[t_, T_] := 
 Log[1/T NIntegrate[(P[t + \[Tau]] - P[\[Tau]])^2, {\[Tau], 0, T - t},
       Method -> "QuasiMonteCarlo"] + 1]/Log[t]

Now we can prepare lists for small, medium and large t:

lst = Table[{t, Kk[t, time]}, {t, 1000, 100000, 3000}];
lst1 = Table[{t, Kk[t, time]}, {t, 100, 10000, 300}]
lst2 = Table[{t, Kk[t, time]}, {t, 10, 100, 5}]

Visualization Figure 1 Therefore at large t function Kk[t] close to constant about 0.75.

  • $\begingroup$ I didn't notice this before. Why have you changed the limit of integration in $M(t)$ from $\tau=0 \to T$ to $\tau=0 \to T-t$ ? It gives the right answer but I can't see why the change is valid. $\endgroup$ Nov 19, 2020 at 7:50
  • $\begingroup$ Actually numerical function $P(t)$ is not defined at $t+\tau>time$. When you try to integrate over this limitation Mathematica automatically use extrapolation. It is why you got answer that K is around 6. Also very important second limitation t<<T. If we check picture I have posted, then we realise that small and large t not valid for this problem. $\endgroup$ Nov 19, 2020 at 10:53

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.