MHB Numerically Solving the Van de Pol Equation to Find Limit Cycle Periods

[Dustinsfl]

s = NDSolve[{x'[t] == \[Mu]*(y[t] - (1/3*x[t]^3 - x[t])), 
    y'[t] == -1/\[Mu]*x[t], x[0] == 8, y[0] == 2}, {x, y}, {t, 150}];

This code will show limit cycles when you parametric plot the Van de Pol Equation. However, I want to find the period T numerically.
I need to add in a piece for mu ranging from 0 - .5 by steps of .05 and 5-50 by steps of 5. After that, I need the code to find the period of the limit cycle. How can I accomplish this?

 

[Ackbach]

To get your $\mu$ stuff going, just enclose your whole procedure inside a for loop. The exact syntax of the Mathematica for loop is

For[start,test,increment,body];

So, for example, you could do this:

MuIncrement = 0.05;
For[Mu=0, Mu <= 0.5, Mu = Mu + MuIncrement,
Solve your ODE problem here, using semi-colons to separate statements
]

MuIncrement = 5;
For[Mu = 5, Mu <= 50, Mu = Mu + MuIncrement,
Solve your ODE problem here, using semi-colons to separate statements
]

As for getting Mathematica to find your limit cycle periods, that's a bit trickier. I would recommend plotting things to see what an approximate solution might be. I'll have to wait on further input until I get on a computer that has Mathematica.

 

[Sudharaka]

s = NDSolve[{x'[t] == \[Mu]*(y[t] - (1/3*x[t]^3 - x[t])), 
    y'[t] == -1/\[Mu]*x[t], x[0] == 8, y[0] == 2}, {x, y}, {t, 150}];

This code will show limit cycles when you parametric plot the Van de Pol Equation. However, I want to find the period T numerically.
I need to add in a piece for mu ranging from 0 - .5 by steps of .05 and 5-50 by steps of 5. After that, I need the code to find the period of the limit cycle. How can I accomplish this?

Hi dwsmith, :)

I posted this question in Google groups and this is the answer that I obtained. I hope you'll find it helpful. :)

Kind Regards,
Sudharaka.

I guess the main problem is the idea to find the period. The rest you may do using, say, the Table function.

Let us look for the period in the case mu=0.5. This is your equation and its solution:

\[Mu] = 0.5;
s = NDSolve[{x'[t] == \[Mu]*(y[t] - (1/3*x[t]^3 - x[t])),
y'[t] == -1/\[Mu]*x[t], x[0] == 8, y[0] == 2}, {x, y}, {t, 0, 150}]{{x -> \!\(\*
TagBox[
RowBox[{"InterpolatingFunction", "[",
RowBox[{
RowBox[{"{",
RowBox[{"{",
RowBox[{"0.`", ",", "150.`"}], "}"}], "}"}], ",", "\<\"<>\"\>"}],
"]"}],
False,
Editable->False]\), y -> \!\(\*
TagBox[
RowBox[{"InterpolatingFunction", "[",
RowBox[{
RowBox[{"{",
RowBox[{"{",
RowBox[{"0.`", ",", "150.`"}], "}"}], "}"}], ",", "\<\"<>\"\>"}],
"]"}],
False,
Editable->False]\)}}

Here we can look at it to make sure the cycle is indeed there:

Plot[Evaluate[x[t] /. s], {t, 0, 20}]

This makes a list of pairs of points {t, x[t]}, in which x[t] has a maximum

lst1 = Last /@
Select[Split[
Table[{t, x[t]} /. s[[1, 1]] /. t -> 0.01*i, {i, 0,
15000}], #1[[2]] < #2[[2]] &], Length[#] > 1 &]

{{6.33, 2.04083}, {12.74, 2.00391}, {19.12, 2.00254}, {25.5,
2.00249}, {31.88, 2.00249}, {38.26, 2.00249}, {44.64,
2.00249}, {51.02, 2.00249}, {57.4, 2.00249}, {63.78,
2.00248}, {70.16, 2.00248}, {76.54, 2.00248}, {82.92,
2.00247}, {89.3, 2.00246}, {95.69, 2.00247}, {102.07,
2.00247}, {108.45, 2.00248}, {114.83, 2.00248}, {121.21,
2.00248}, {127.59, 2.00249}, {133.97, 2.00249}, {140.35,
2.00249}, {146.73, 2.00249}, {150., -1.99652}}

This calculates differences between the times corresponding to successive points of maximum:

Lst2=Differences[Transpose[lst1][[1]]]

{6.41, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, \
6.38, 6.38, 6.39, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, 6.38, \
3.27}

The few first values should be dropped, since the system may have not fallen on the cycle yet, while the very last should be omitted, since you may have interrupted solving your system in an arbitrary moment of time. And this is your period:

Drop[Drop[lst2, 3], -1] // Mean6.38053

Few notes:

1) The precision in the obtained period is defined by what we choose in the Table operator, namely, t -> 0.01*i yields you the precision of 0.01 (and not 0.00001 as may be erroneously concluded from the final result). But you can increase precision, if needed.

2) I do not see, a good way of how to automate the search of the number of items to drop. Something like this, may be?

myRound[x_] := Round[100.*x]/100.;
lst3 = Split[lst2, myRound[#1] == myRound[#2] &]

{{6.376}, {6.313}, {6.297}, {6.293, 6.289, 6.289, 6.287, 6.288, 6.287,
6.287, 6.288, 6.287, 6.287, 6.287, 6.287, 6.287, 6.287, 6.287,
6.288, 6.287, 6.287, 6.287}, {5.44}}

Select[lst3, Length[#] > 1 &] // Flatten // Mean

6.28768

Here, however, I had to set t -> 0.001*i in the Table, otherwise it was too rough. So one needs to check it once by eye. I hope you will solve this question.