Re: Ordinary Differential Equation: Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg128397] Re: Ordinary Differential Equation: Mathematica
- From: Alexei Boulbitch <Alexei.Boulbitch at iee.lu>
- Date: Tue, 16 Oct 2012 03:20:55 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-newout@smc.vnet.net
- Delivered-to: mathgroup-newsend@smc.vnet.net
Hi Everyone, :)
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 above 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 to 0.5 by steps of 0.05 and 5 to 50 by steps of 5. After that, I need the code to find the period of the limit cycle. How can I accomplish this?
Kind Regards,
Sudharaka.
Hi, 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] // Mean
6.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.
Have fun, Alexei
Alexei BOULBITCH, Dr., habil.
IEE S.A.
ZAE Weiergewan,
11, rue Edmond Reuter,
L-5326 Contern, LUXEMBOURG
Office phone : +352-2454-2566
Office fax: +352-2454-3566
mobile phone: +49 151 52 40 66 44
e-mail: alexei.boulbitch at iee.lu