# Not able to reproduce the result of this paper to test if a system is chaotic or nonchaotic

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

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$$

And, finally using $$p(t)$$ calculate

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.

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}]]/
Log[t]
K[Exp[10], Exp[10]]
• 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. Commented Nov 17, 2020 at 11:23
• @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? Commented Nov 17, 2020 at 12:55
• @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. Commented Nov 17, 2020 at 13:42
• 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. Commented Nov 17, 2020 at 13:46
• Now use second restriction underlined in the paper that $t<<T$ and use Log[M[t]+1]/Log[t] as they recommended. Commented Nov 17, 2020 at 13:51

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 Therefore at large t function Kk[t] close to constant about 0.75.

• 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. Commented Nov 19, 2020 at 7:50
• 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. Commented Nov 19, 2020 at 10:53