(* by Nasser M. Abbasi, version: April 26, 2013 *)
Manipulate[
(
(*If Q, the quality factor is zero,
then we assume no damping will be used*)
(*This way there is no need to add an extra dialog asking if \
damping is to be used or not*)
isDamped = If[Abs@q <= $MachineEpsilon, False, True];
(* start of the finite state machine,
where all the logic of the program in included. 5 states,
8 events *)
(* events are set using the second argument of dynamics right at \
the control level below *)
(* this is a way to simulate event driven UI with callbacks *)
Which[
(*------- STATE RESET ---------------*)
state == "reset",
(
currentTime = 0;
tick = 0;
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, x1,
x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime];
If[bifurcationPlot === {},
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress]
];
state = "pause"
),
(*------- STATE PAUSE ---------------*)
state == "pause",
(
Which[
lastEvent == "no_event", state = "pause",
lastEvent == "pause_button", (lastEvent = "no_event";
state = "pause"),
lastEvent == "reset_button", (lastEvent = "no_event";
state = "reset"; tick += DEL),
lastEvent ==
"initial_conditions_changed", (lastEvent = "no_event";
state = "reset"; tick += DEL),
lastEvent == "time_series_choice", (lastEvent = "no_event";
bottomPlot = If[plotChoice == 1,
makeTimeSeriesPlot[phaseData, plotImageSize, nSteps,
currentTime],
makePowerSpectrumPlot[phaseData, plotImageSize, nSteps, omega,
spectrumScale, currentTime]
]),
lastEvent == "step_button",
(
lastEvent = "no_event";
{currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3} =
makeOneStep[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, nSteps, dt, phasePlot, bottomPlot,
track, maxRunTime, x1, x2, x3]
),
lastEvent == "run_button", (lastEvent = "no_event";
state = "running"; tick += DEL),
lastEvent == "fast_button", (lastEvent = "no_event";
state = "fast mode"; tick += DEL),
lastEvent == "update_bifurcation_button",
(lastEvent = "no_event";
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress]
),
lastEvent == "test_case_changed",
(
lastEvent = "no_event";
{torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals,
x10, isDamped, x20} = setParameters[testCase];
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress];
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime]
)
]
),
(*------- STATE RUNNING ---------------*)
state == "running",
(
If[poincareMap === {},
(
currentTime = 0;
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime]
)
];
If[bifurcationPlot === 0,
(
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress]
)
];
Which[
lastEvent == "initial_conditions_changed",
(
lastEvent = "no_event";
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime];
tick += DEL
),
lastEvent == "no_event" || lastEvent == "time_series_choice",
(
If[lastEvent == "time_series_choice", lastEvent = "no_event"];
If[currentTime + dt <= maxRunTime,
(
{currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3} =
makeOneStep[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, nSteps, dt, phasePlot,
bottomPlot, track, maxRunTime, x1, x2, x3];
tick += DEL
),
(
state = "pause"
)
]
),
lastEvent == "step_button",
(
lastEvent = "no_event";
If[currentTime + dt <= maxRunTime,
(
{currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3} =
makeOneStep[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, nSteps, dt, phasePlot,
bottomPlot, track, maxRunTime, x1, x2, x3];
state = "pause"
)
]
),
lastEvent == "pause_button", (lastEvent = "no_event";
state = "pause"),
lastEvent == "fast_button", (lastEvent = "no_event";
state = "fast mode" ; tick += DEL),
lastEvent == "reset_button", (lastEvent = "no_event";
state = "reset"; tick += DEL),
lastEvent == "update_bifurcation_button",
(
lastEvent = "no_event";
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress];
tick += DEL
),
lastEvent == "test_case_changed",
(
lastEvent = "no_event";
{torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals,
x10, isDamped, x20} = setParameters[testCase];
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress];
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime];
tick += DEL
)
]
),
(*------- STATE FAST_MODE ---------------*)
state == "fast mode",
(
If[poincareMap === {},
(
currentTime = 0;
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime]
)
];
If[bifurcationPlot === {},
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress]
];
Which[
lastEvent == "no_event" ||
lastEvent == "initial_conditions_changed" ||
lastEvent == "time_series_choice" || lastEvent == "fast_button",
(
If[Not[lastEvent == "no_event"], lastEvent = "no_event"];
currentTime = maxRunTime;
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, currentRunMaxSpeed} =
initializeSystem[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel, plotChoice, phaseData,
spectrumScale, currentTime, dt, maxRunTime]
),
lastEvent == "pause_button", (lastEvent = "no_event";
state = "pause"),
lastEvent == "reset_button", (lastEvent = "no_event";
state = "reset"; tick += DEL),
lastEvent == "run_button", (lastEvent = "no_event"),
lastEvent == "step_button", (lastEvent = "no_event"),
lastEvent == "update_bifurcation_button",
(lastEvent = "no_event";
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress]
),
lastEvent == "test_case_changed",
(
lastEvent = "no_event";
{torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals,
x10, isDamped, x20} = setParameters[testCase];
bifurcationPlot =
makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase,
x10, x20, isDamped, plotImageSize,
Unevaluated@bifurcationProgress];
tick += DEL
)
]
)
];
generateFinalDisplay[poincareMap, phasePlot, bifurcationPlot, track,
currentTime, plotImageSize, x1, x2, x3, bottomPlot,
currentRunMaxSpeed, torqueAmplitude]
),
(*--- START OF CONTROLS ---*)
Grid[{
{Grid[{(* BLOCK 1 *)
{Grid[{
{Button[
Text[Style["play", 14]], (lastEvent = "run_button";
tick += DEL), ImageSize -> {75, 30}],
Button[Text[
Style["pause", 14]], (lastEvent = "pause_button";
tick += DEL), ImageSize -> {75, 30}],
Button[Text[Style["step", 14]], (lastEvent = "step_button";
tick += DEL), ImageSize -> {75, 30}]}
}, Spacings -> {.2, .1}, Alignment -> Left]
},
{Grid[{
{Button[
Text[Style["reset", 14]], (lastEvent = "reset_button";
tick += DEL), ImageSize -> {75, 30}],
Button[Text[Style["fast", 14]], (lastEvent = "fast_button";
tick += DEL), ImageSize -> {75, 30}],
Framed[Dynamic@
Graphics[Text@Style[state, 12, Italic],
ImageSize -> {75, 30}],
FrameStyle -> Directive[Thickness[.005], Gray],
FrameMargins -> -1, BaselinePosition -> Baseline,
ImageMargins -> 0,
ContentPadding -> True]} (*the state of the system*)
}, Spacings -> {.2, .1}, Alignment -> Left]}
}, Spacings -> {.2, .1}, Alignment -> Center
]
},
{Grid[{(* BLOCK 2 *)
{Text@Style["maximum duration", 11], Spacer[4],
Manipulator[
Dynamic[maxRunTime, (maxRunTime = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {1,
500, 1}, ImageSize -> Tiny, ContinuousAction -> False],
Spacer[4], Dynamic@Text@Style[padIt3[maxRunTime, {3, 0}], 11]
},
{Text@
Style[TraditionalForm[HoldForm[\[CapitalDelta]\[Tau]]], 11],
Spacer[4],
Manipulator[
Dynamic[dt, (dt = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {0.01,
1, 0.01}, ImageSize -> Tiny, ContinuousAction -> False],
Spacer[4], Dynamic@Text@Style[padIt2[dt, {3, 3}], 11],
SpanFromLeft}
}, Frame -> None, Spacings -> {0, 0.1}, Alignment -> Center
]
},
{Grid[{(* BLOCK 3 *)
{Row[{Style["\[Theta]''(\[Tau]) + \[Theta]'(\[Tau])/", "TR",
12], Style["Q", Italic, "TR", 12],
Style["+ sin \[Theta](\[Tau])) \[Equal] ", "TR", 12] ,
Style["a", Italic, "TR", 12],
Style[" cos \[Phi](\[Tau])", "TR", 12]}]}
}]
},
{Grid[{(* BLOCK 4 *)
{Style["Q", Italic, "TR"], Spacer[2],
Manipulator[
Dynamic[q, (q = #; {lastEvent = "initial_conditions_changed",
tick += DEL}; #) &], {0, 4, 0.01}, ImageSize -> Tiny,
ContinuousAction -> False], Spacer[2],
Dynamic@Text@Style[padIt2[q, {3, 3}], 12], "", ""},
{ Text@Style["\[Omega] = d\[Phi]/d\[Tau]", 12], Spacer[2],
Manipulator[
Dynamic[omega, (omega = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {-4,
4, 0.01}, ImageSize -> Tiny, ContinuousAction -> False],
Spacer[2], Dynamic@Text@Style[padIt1[omega, {3, 3}], 11], "",
""},
{Text@Style["a", Italic, 12], Spacer[2],
Manipulator[
Dynamic[torqueAmplitude, (torqueAmplitude = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {0, 2,
0.01}, ImageSize -> Tiny, ContinuousAction -> False],
Spacer[2],
Dynamic@Text@Style[padIt2[torqueAmplitude, {3, 3}], 11], "", ""}
}, Frame -> None, Spacings -> {0.2, 0.6}, Alignment -> Center]
},
{Grid[{(* BLOCK 5 *)
{Item[Text@Style["initial conditions", 12],
Alignment -> Center], SpanFromLeft},
{ Text@Style["\[Theta] (0)", 12], Spacer[2],
Manipulator[
Dynamic[x10, (x10 = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {-N@
Pi, N@Pi, 0.01}, ImageSize -> Tiny,
ContinuousAction -> False], Spacer[2],
Dynamic@Text@Style[padIt1[x10, {6, 4}], 11], Spacer[1],
Text@Style["radian", 10]},
{Text@Style["d\[Theta]/d\[Tau] (0)", 12], Spacer[2],
Manipulator[
Dynamic[x20, (x20 = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {-Pi,
Pi, 0.01}, ImageSize -> Tiny, ContinuousAction -> False],
Spacer[2], Dynamic@Text@Style[padIt1[x20, {6, 4}], 11], "", ""},
{Text@Style["\[Phi] (0)", 12], Spacer[2],
Manipulator[
Dynamic[phase, (phase = #; {lastEvent =
"initial_conditions_changed", tick += DEL}; #) &], {-N@
Pi, N@Pi, 0.01}, ImageSize -> Tiny,
ContinuousAction -> False], Spacer[2],
Dynamic@Text@Style[padIt1[phase, {5, 3}], 11], Spacer[1],
Text@Style["radian", 11]}
}, Spacings -> {.4, 0.2}, Alignment -> Left
]},
{Grid[{(* BLOCK 6 *)
{Item[
Text@Row[{Style["select range of parameter ", 12],
TraditionalForm[HoldForm[a]]}], Alignment -> Center],
SpanFromLeft}, {Item[Text@Style["from", 11],
Alignment -> Left],
Control[{{aStart, 0.9, ""}, 0, 2, 0.01, ImageSize -> Tiny,
ImagePadding -> 0}],
Dynamic@Text@Style[padIt2[aStart, {3, 3}], 11]},
{Item[Text@Style["length", 11], Alignment -> Left],
Control[{{aLen, 0.6, ""}, 0.01, 2, 0.01, ImageSize -> Tiny,
ImagePadding -> 0}],
Dynamic@Text@Style[padIt2[aLen, {3, 3}], 11]},
{Item[Text@Style["intervals", 11], Alignment -> Left],
Control[{{aIntervals, 20, ""}, 5, 300, 1, ImageSize -> Tiny,
ImagePadding -> 0}],
Dynamic@Text@Style[padIt3[Rationalize@aIntervals, {4, 0}], 11]},
{Button[
Text[Style["generate", 11]], {lastEvent =
"update_bifurcation_button"; tick += DEL;
Method -> "Queued"}, ImageSize -> {65, 25},
BaselinePosition -> Bottom],
Dynamic@ProgressIndicator[bifurcationProgress, {0, 1},
ImageSize -> {105, 23}, ImageMargins -> 0,
BaselinePosition -> Bottom],
Dynamic@Graphics[
Text@Row[{
Style[padIt2[bifurcationProgress*100, {4, 2}] , 11],
Style["%", 11]}], ImageSize -> {40, 20},
BaselinePosition -> Bottom]}
}, Frame -> None, Spacings -> {.4, 0.1}, Alignment -> Center]
},
{Grid[{(* BLOCK 8 *)
{Text@Style["preconfigured test cases", 10]}, {PopupMenu[
Dynamic[
testCase, (testCase = #; {lastEvent = "test_case_changed" ,
tick += DEL}; #) &],
{
"double periodic #1" ->
Style["double periodic, on edge of chaotic", 10],
"periodic, A=0.9,Q=2,omega=2/3" ->
Style["periodic, A=0.9,Q=2,omega=2/3", 10],
"periodic, A=1.35,Q=2,omega=2/3" ->
Style["periodic, A=0.9,Q=2,omega=2/3", 10],
"period doubling, A=1.12,Q=2,omega=2/3" ->
Style["period doubling, A=1.12,Q=2,omega=2/3", 10],
"period doubling, A=1.07,Q=2,omega=2/3" ->
Style["period doubling, A=1.07,Q=2,omega=2/3", 10],
"period doubling, A=1.45,Q=2,omega=2/3" ->
Style["period doubling, A=1.45,Q=2,omega=2/3", 10],
"period quadrupling, A=1.47,Q=2,omega=2/3" ->
Style["period quadrupling, A=1.47,Q=2,omega=2/3", 10],
"chaotic, A=1.15,Q=2,omega=2/3" ->
Style["chaotic, A=1.15,Q=2,omega=2/3", 10],
"chaotic, A=1.5,Q=2,omega=2/3" ->
Style["chaotic, A=1.5,Q=2,omega=2/3", 10],
"chaotic, A=1.5,Q=2,omega=2/3,phase=3.14" ->
Style["chaotic, A=1.5,Q=2,omega=2/3,phase=3.14", 10],
"chaotic, Bifurcation #1" ->
Style["chaotic, Bifurcation #1", 10],
"chaotic, Bifurcation #2" ->
Style["chaotic, Bifurcation #2", 10],
"chaotic, Bifurcation #3" ->
Style["chaotic, Bifurcation #3", 10]}, ImageSize -> All,
ContinuousAction -> False]
}
}
]
},
{Grid[{(* BLOCK 9 *)
{Text@Style["plot type ", 10], Spacer[5],
RadioButtonBar[
Dynamic[plotChoice, (plotChoice = #; {lastEvent =
"time_series_choice", tick += DEL}; #) &], {1 ->
Text@Style["time series", 11],
2 -> Text@Style["power spectrum", 11]}]}, {Text@
Style["spectrum scale", 10], Spacer[5],
RadioButtonBar[
Dynamic[spectrumScale, (spectrumScale = #; {lastEvent =
"time_series_choice", tick += DEL}; #) &], {1 ->
Text@Style["standard ", 11], 2 -> Text@Style["log", 11]},
Enabled -> Dynamic[plotChoice == 2]]}
}, Frame -> None, Spacings -> {0.1, 0.2}, Alignment -> Left]
}
}, Frame -> All, FrameStyle -> Directive[Thickness[.005], Gray],
Spacings -> {0.2, 0.6}, Alignment -> Center
],
{{DEL, $MachineEpsilon}, None},
{{currentRunMaxSpeed, 0}, None},
{{tick, 0}, None},
{{currentTime, 300}, None},
{{plotChoice, 1}, None},
{{spectrumScale, 1}, None},
{{bifurcationProgress, 0}, None},
{{omega, 0.6667}, None},
{{q, 2}, None},
{{torqueAmplitude, 1.5}, None},
{{phase, 0}, None},
{{x10, 0.6184}, None},
{{x20, 0}, None},
{{dt, 0.15}, None},
{{x1, 0}, None},
{{x2, 0}, None},
{{x3, 0}, None},
{{maxRunTime, 300}, None},
{{isDamped, True}, None},
{{testCase, "periodic, A=0.9,Q=2,omega=2/3"}, None},
{{state, "fast mode"}, None},
{{lastEvent, "no_event"}, None},
{{track, {Blue, PointSize[0.03], Point[ {0, 0}]}}, None},
{{phaseData , Table[{0, 0, 0}, {500*101}]},
None}, (*DO NOT CHANGE, DEPENDS ON FIXED SETTINGS *)
{{poincareMap, {}}, None},
{{phasePlot, {}}, None},
{{bottomPlot, {}}, None},
{{nSteps, 0}, None},
{{plotImageSize, {340, 340}}, None},
{{axesLabel, {Text@Style["\[Theta]", 9],
Text@Style[
Row[{" ", Style["d", Italic], "\[Theta]/", Style["d", Italic],
"\[Tau]"}], 9]}}, None},
{{bifurcationPlot, {}}, None},
TrackedSymbols -> {tick}, (*has to be -> and NOT :> else won't work *)
SynchronousUpdating -> False,
ContinuousAction -> False,
SynchronousInitialization -> False,
ControlPlacement -> Left,
Alignment -> Center, ImageMargins -> 1, FrameMargins -> 1,
Initialization :>
(
(*----------------------------------------------------------------------*)
padIt1[v_?(NumericQ[#] && Im[#] == 0 &), f_List] :=
AccountingForm[Chop[N@v] , f, NumberSigns -> {"-", "+"},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*-----------------------------------------------------------------------*)
padIt2[v_?(NumericQ[#] && Im[#] == 0 &), f_List] :=
AccountingForm[Chop[N@v] , f, NumberSigns -> {"", ""},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*--------------------------------------------------------------------*)
padIt3[v_?(NumericQ[#] && Im[#] == 0 &), f_List] :=
AccountingForm[v , f, NumberSigns -> {"", ""},
NumberPadding -> {"0", "0"}, SignPadding -> True,
NumberPoint -> ""];
(*------------------------------------------------------------------------*)
getBobPosition[len_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
angle_?(NumericQ[#] && Im[#] == 0 &)] := Module[{bob},
bob = {len Sin[angle], -len Cos[angle]};
{ {Red, Thick, Style[Line[{{0, 0}, bob}], Antialiasing -> True]},
{Blue, Opacity[1], PointSize[0.06],
Style[Point[bob], Antialiasing -> True]} }
];
(*-----------------------------------------------------------------------*)
(* keep angle from -180...+180 to make phase and poincare plots easier to make*)
normalize[angle_?(NumericQ[#] && Im[#] == 0 &)] :=
angle - 2 Pi Floor[(angle + Pi)/(2 Pi)];
(*---------------------------------------------------------------------------------*)
solve[q_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*quality factor*)
a_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*driver amplitude*)
omega_?(NumericQ[#] &&
Im[#] == 0 &), (*driver frequency*)
phase_?(NumericQ[#] && Im[#] == 0 &), (*driver phase*)
from_?(NumericQ[#] &&
Im[#] ==
0 &), (*time to start integrating the ode*)
to_?(NumericQ[#] && Im[#] == 0 &), (*end time *)
x10_?(NumericQ[#] &&
Im[#] ==
0 &), (*initial condition for angle at time=0*)
x20_?(NumericQ[#] &&
Im[#] ==
0 &), (*initial condition for speed at time=0*)
isDamped_?(Element[#, Booleans] &), (*is Q used?*)
accuracyGoal_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
precisionGoal_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &)] :=
Module[ {ic, eq1, eq2, eq3, x1, x2, x3, t, sol, ndsolveOptions},
ndsolveOptions = {MaxSteps -> Infinity,
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}},
AccuracyGoal -> accuracyGoal, PrecisionGoal -> precisionGoal};
eq1 = x1'[t] == x2[t];
If[isDamped,
eq2 = x2'[t] == -1/q x2[t] - Sin[x1[t]] + a Cos[x3[t]],
eq2 = x2'[t] == -Sin[x1[t]] + a Cos[x3[t]]
];
eq3 = x3'[t] == omega;
ic = {x1[0] == x10, x2[0] == x20, x3[0] == phase};
sol =
First@NDSolve[
Flatten@{eq1, eq2, eq3, ic}, {x1, x2, x3}, {t, from, to},
Sequence@ndsolveOptions];
{x1 /. sol, x2 /. sol, x3 /. sol}
];
(*---------------------------------------------------------------------------------*)
makePoincreSection[
q_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*quality factor*)
torqueAmplitude_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*driver amplitude*)
omega_?(NumericQ[#] &&
Im[#] == 0 &), (*driver frequency*)
phase_?(NumericQ[#] &&
Im[#] == 0 &), (*driver phase*)
x10_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for angle at time=0*)
x20_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for speed at time=0*)
isDamped_?(Element[#, Booleans] &), (*is Q used?*)
plotImageSize_List, axesLabel_List] :=
Module[{x1, x2, x3, i, data, accuracyGoal = 4, precisionGoal = 4,
pointSize = 0.01, driverPeriod, initialDriverPeriod = 150},
If[Abs@omega <= $MachineEpsilon ||
torqueAmplitude <= $MachineEpsilon,
(
pointSize = 0.05;
data = {{x10, x20}}
),
(
driverPeriod = 2. Pi/Abs@omega;
{x1, x2, x3} =
solve[q, torqueAmplitude, omega, phase,
initialDriverPeriod*driverPeriod, 1000*driverPeriod, x10,
x20, isDamped, accuracyGoal, precisionGoal];
data =
Table[{normalize[x1[i]], x2[i]}, {i,
initialDriverPeriod*driverPeriod, 1000*driverPeriod,
driverPeriod}]
)
];
ListPlot[data, AspectRatio -> 1,
Ticks -> {{-Pi, {-Pi/2, "-\[Pi]/2"}, 0, {Pi/2, "\[Pi]/2"}, Pi},
Automatic},
PlotStyle -> {Blue, Directive[PointSize[pointSize]]},
PerformanceGoal -> "Quality", ImageSize -> plotImageSize/2,
PlotLabel -> Text@Style["Poincaré section", 12, Bold],
PerformanceGoal -> "Quality", TicksStyle -> Directive[7],
AxesLabel -> {Text@Style["\[Theta]", 9],
Text@Style[
Row[{" sampled ", Style["d", Italic], "\[Theta]/",
Style["d", Italic], "\[Tau]"}], 9]}
]
];
makePhasePlot[data_List,(*phase data*)
plotImageSize_List, axesLabel_List,
nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &)] := Module[{},
ListPlot[data[[1 ;; nStep, 2 ;; 3]], Joined -> False,
PlotStyle -> {Red, Directive[PointSize[0.003]]},
PerformanceGoal -> "Quality",
PlotLabel -> Text@Style["phase portrait", 12, Bold],
AxesLabel -> axesLabel, ImageSize -> plotImageSize/2,
AspectRatio -> 1,
PlotRange -> {{-1.1 Pi, 1.1 Pi}, {-1.2 Pi, 1.2 Pi}},
Ticks -> {{-Pi, {-Pi/2, "-\[Pi]/2"}, 0, {Pi/2, "\[Pi]/2"}, Pi},
Automatic}, PlotStyle -> {Blue, Directive[PointSize[0.003]]},
TicksStyle -> Directive[8],
PerformanceGoal -> "Quality"
]
];
(*---------------------------------------------------------------------------------*)
generateFinalDisplay[poincareMap_Graphics,
phasePlot_Graphics,
bifurcationPlot_Graphics,
track_List,
currentTime_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),
plotImageSize_List, x1_, x2_, x3_,
bottomPlot_Graphics,
currentRunMaxSpeed_?(NumericQ[#] && Im[#] == 0 &),
torqueAmplitude_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &) (*driver amplitude*)] :=
Module[{g, options, arrow, torqueArrow, speed = x2[currentTime],
angle = x1[currentTime], angleOffset, torqueAngleOffset,
gTorque},
angleOffset = (Abs@speed/currentRunMaxSpeed)*0.35*Pi;
If[angleOffset > $MachineEpsilon,
arrow = If[speed < 0,
arcArrow[{0, 0}, 1, angle - Pi/2 - angleOffset, angle - Pi/2 ,
True], arcArrow[{0, 0}, 1,
angle - Pi/2, (angle + angleOffset) - Pi/2 , False]
]
];
If[torqueAmplitude <= $MachineEpsilon,
gTorque = {},
(
torqueAngleOffset = Abs@Cos[x3[currentTime]]*Pi;
If[torqueAngleOffset > 0,
(
torqueArrow = If[Cos[x3[currentTime]] < 0,
arcArrow[{0, 0}, 0.25, Pi - torqueAngleOffset , Pi, True],
arcArrow[{0, 0}, 0.25, -Pi, -Pi + torqueAngleOffset , False]
];
gTorque = {Opacity[1], Thickness[0.010], Arrowheads[.08],
torqueArrow}
)
];
)
];
g = Graphics[{
{Dotted, RGBColor[65/255, 105/255, 225/255], Thickness[0.01],
Circle[{0, 0}, 1]},
{Red, PointSize[0.03], Point[{0, 0}]},
getBobPosition[1, angle],
If[
angleOffset > $MachineEpsilon, {Thick,
RGBColor[205/255, 92/255, 92/255], Arrowheads[.07],
arrow}, {}],
gTorque
}];
options = {ImageMargins -> 1, Axes -> False,
PlotRange -> {{-1.2, 1.2 }, {-1.2 , 1.2}}, AspectRatio -> 1,
Background -> White, TicksStyle -> Directive[8],
ImagePadding -> 4, ImageSize -> plotImageSize/2,
AxesStyle -> Directive[{Blue, Thickness[0.01]}]};
Grid[{
{poincareMap, Show[phasePlot, Graphics[track]]},
{bifurcationPlot, Show[g, Sequence@options]},
{bottomPlot, SpanFromLeft}
}, Frame -> All,
FrameStyle -> Directive[Thickness[.005], Gray],
Spacings -> {0.5, 0.4}
]
];
(*---------------------------------------------------------------------------------*)
makeBifurcationPlot[
aStart_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*driver amplitude starting value*)
aLen_?(NumericQ[#] && Im[#] == 0 &&
Re[#] > 0 &), (*driver amplitude range of values length*)
aIntervals_?(IntegerQ[#] && Im[#] == 0 &&
Re[#] > 0 &),(*how many intervals to do *)
q_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*quality factor*)
omega_?(NumericQ[#] &&
Im[#] == 0 &), (*driver frequency*)
phase_?(NumericQ[#] &&
Im[#] == 0 &), (*driver phase*)
x10_?(NumericQ[#] &&
Im[#] ==
0 &), (*initial condition for angle at time=
0*)x20_?(NumericQ[#] &&
Im[#] ==
0 &), (*initial condition for speed at time=
0*)
isDamped_?(Element[#, Booleans] &),
plotImageSize_List,
bifurcationProgress_ (*by reference*)] :=
Module[{driverPeriod, timeValues, valuesOfDriverAmplitude, x1, x2,
x3, i, j, p, data, initialDriverPeriod = 150,
finalDriverPeriod = 250, plotOptions, total, count,
currentTorqueAmplitude, pointSize = 0.004},
plotOptions = {Joined -> False,
PlotStyle -> {Red, Directive[PointSize[pointSize]]},
PerformanceGoal -> "Quality",
PlotLabel -> Text@Style["bifurcation map", 12, Bold],
ImageSize -> plotImageSize/2, AspectRatio -> 0.85,
PlotRange -> {{aStart, aStart + aLen}, All},
PlotStyle -> {Blue, Directive[PointSize[0.003]]},
TicksStyle -> Directive[8],
AxesLabel -> {Text@Style["a", Italic, 9],
Text@Style[
Row[{" sampled ", Style["d", Italic], "\[Theta]/",
Style["d", Italic], "\[Tau]"}], 9]},
PerformanceGoal -> "Quality"};
If[Abs@omega <= $MachineEpsilon,
(
data = {{0, 0}};
p = ListPlot[data, Sequence@plotOptions]
)
,
(
driverPeriod = 2. Pi/Abs@omega;
valuesOfDriverAmplitude =
Range[aStart, aStart + aLen, aLen/aIntervals];
timeValues =
Range[initialDriverPeriod*driverPeriod,
finalDriverPeriod*driverPeriod, driverPeriod ];
data =
Table[0, {Length[valuesOfDriverAmplitude]}, {Length[
timeValues ]}];
total = Length[timeValues]*Length[valuesOfDriverAmplitude];
bifurcationProgress = 0.0;
count = 0.0;
Do[
(
{x1, x2, x3} =
solve[q, valuesOfDriverAmplitude [[i]], omega, phase,
initialDriverPeriod*driverPeriod, (finalDriverPeriod)*
driverPeriod, x10, x20, isDamped, 4, 4];
currentTorqueAmplitude = valuesOfDriverAmplitude[[i]];
Do[
(
data[[i, j]] = {currentTorqueAmplitude,
x2[timeValues[[j]] ]};
count += 1;
bifurcationProgress = count/total
), {j, 1, Length[timeValues]}
]
), {i, 1, Length[valuesOfDriverAmplitude]}
];
p = ListPlot[Flatten[data, 1], Sequence@plotOptions]
)
];
bifurcationProgress = 0;
p
];
(*---------------------------------------------------------------------------------*)
makeTimeSeriesPlot[data_List, plotImageSize_List,
nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
currentTime_] := Module[{},
ListPlot[
Transpose[{data[[1 ;; nStep, 1]], data[[1 ;; nStep, 3]]}],
Joined -> True, GridLines -> Automatic,
PlotStyle -> {Red, Directive[PointSize[0.001]]},
PerformanceGoal -> "Quality",
ImageSize -> {370, .44*plotImageSize}, AspectRatio -> .30,
PlotStyle -> {Blue, Directive[PointSize[0.003]]},
AxesOrigin -> {0, 0}, ImagePadding -> {{30, 2}, {15, 35}},
PlotRange -> {{0, Automatic}, All},
PerformanceGoal -> "Quality", Frame -> True,
FrameTicksStyle -> Directive[7], Axes -> None,
FrameLabel -> {{None, None}, {None,
Text@Row[{Style["time series", 11, Bold], Spacer[5],
Style[Row[{Style["d", Italic], "\[Theta]/",
Style["d", Italic],
"\[Tau] vs. \[Tau] (dimensionless)"}], 9], Spacer[5],
Style["time: ", 10],
Text@Style[padIt2[currentTime, {5, 2}], 11]}]}}
]
];
(*---------------------------------------------------------------------------------*)
makePowerSpectrumPlot[data_List, plotImageSize_List,
nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
omega_?(NumericQ[#] && Im[#] == 0 &),
spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
currentTime_] := Module[{list, len},
list =
Fourier[data[[1 ;; nStep, 3]], FourierParameters -> {1, -1}];
list = If[Abs[omega] > 2*$MachineEpsilon,
(Abs@list)^2/(nStep*Abs@omega),
(Abs@list)^2/nStep
];
len = Length[list];
If[len > 1, list = list[[1 ;; Floor[Length[list]/2] ]]];
ListPlot[If[spectrumScale == 1, list, Log[10, list]],
Joined -> True, GridLines -> Automatic,
PlotStyle -> {Red, Directive[PointSize[0.001]]},
PerformanceGoal -> "Quality",
ImageSize -> {370, 0.44*plotImageSize}, AspectRatio -> .30,
PlotStyle -> {Blue, Directive[PointSize[0.003]]},
AxesOrigin -> {0, 0}, PlotRange -> {{0, Automatic}, All},
ImagePadding -> {{30, 2}, {15, 35}},
PerformanceGoal -> "Quality", Frame -> True,
FrameTicksStyle -> Directive[8], Axes -> None,
RotateLabel -> False,
FrameLabel -> {{None, None}, {None,
Text@Row[{Style["power spectrum of angular velocity", 11,
Bold], Spacer[5], Style[" time: ", 10],
Text@Style[padIt2[currentTime, {5, 2}], 11]}]} }
]
];
(*---------------------------------------------------------------------------------*)
setParameters[testCase_String] :=
Module[{a, omega, phase, aStart, aLen, aIntervals, x10, isDamped,
x20, q},
Which[
testCase == "periodic, A=0.9,Q=2,omega=2/3", (a = 0.9; q = 2;
omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6;
aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase == "periodic, A=1.35,Q=2,omega=2/3", (a = 1.35; q = 2;
omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.7;
aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase == "period doubling, A=1.12,Q=2,omega=2/3", (a = 1.12;
q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.7;
aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase == "period doubling, A=1.07,Q=2,omega=2/3", (a = 1.07;
q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6;
aIntervals = 100; x10 = 45*Pi/180; isDamped = True; x20 = 0),
testCase == "period doubling, A=1.45,Q=2,omega=2/3", (a = 1.45;
q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6;
aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase ==
"period quadrupling, A=1.47,Q=2,omega=2/3", (a = 1.47; q = 2;
omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6;
aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase == "chaotic, A=1.15,Q=2,omega=2/3",
(a = 1.15; q = 2; omega = 2/3; phase = 0; aStart = 0.9;
aLen = 0.7; aIntervals = 50; x10 = 45*Pi/180.; isDamped = True;
x20 = 0),
testCase == "chaotic, A=1.5,Q=2,omega=2/3",
(a = 1.5; q = 2; omega = 2/3; phase = 0; aStart = 0.9;
aLen = 0.7; aIntervals = 50; x10 = 45*Pi/180.; isDamped = True;
x20 = 0),
testCase == "chaotic, A=1.5,Q=2,omega=2/3,phase=3.14", (a = 1.5;
q = 2; omega = 2/3; phase = 3.14; aStart = 0.9; aLen = 0.7;
aIntervals = 50; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
testCase === "chaotic, Bifurcation #1",
(a = 1.5; q = 2; omega = 2/3; phase = 0; aStart = 0.9;
aLen = 0.2; aIntervals = 70; x10 = 45*Pi/180.; isDamped = True;
x20 = 0),
testCase === "chaotic, Bifurcation #2",
(a = 1.15; q = 4; omega = 2/3; phase = 0; aStart = 0.9;
aLen = 0.8; aIntervals = 80; x10 = 45*Pi/180.; isDamped = True;
x20 = 0),
testCase === "chaotic, Bifurcation #3",
(a = 1.5; q = 2; omega = -1.33; phase = 1.8; aStart = 0.9;
aLen = 0.6; aIntervals = 150; x10 = 50*Pi/180.;
isDamped = True; x20 = 1.642),
testCase === "double periodic #1", (a = 1.08; q = 2;
omega = N[2/3]; phase = 0; aStart = 1.050; aLen = 0.08;
aIntervals = 300; x10 = 0; isDamped = True; x20 = 0)
];
{a, q, omega, phase, aStart, aLen, aIntervals, x10, isDamped, x20}
];
(*---------------------------------------------------------------------------------*)
initializeSystem[
q_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),(*quality factor*)
torqueAmplitude_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &), (*driver amplitude*)
omega_?(NumericQ[#] &&
Im[#] == 0 &), (*driver frequency*)
phase_?(NumericQ[#] &&
Im[#] == 0 &), (*driver phase*)
x10_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for angle at time=0*)
x20_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for speed at time=0*)
isDamped_?(Element[#, Booleans] &),
plotImageSize_List,
axesLabel_List,
plotChoice_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
pphaseData_,
spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
currentTime_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),
dt_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
maxRunTime_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &)] :=
Module[{phaseData = pphaseData, nSteps, poincareMap, phasePlot,
bottomPlot, track, x1, x2, x3, i, range, maxSpeed},
(*solve first to find maximum speed over the whole run,
needed for plotting relative arrow*)
{x1, x2, x3} =
solve[q, torqueAmplitude, omega, phase, 0, maxRunTime, x10, x20,
isDamped, 4, 4];
maxSpeed = Max@Abs@Table[x2[i], {i, 0, maxRunTime, dt}];
{x1, x2, x3} =
solve[q, torqueAmplitude, omega, phase, 0, currentTime + 1, x10,
x20, isDamped, 4, 4];
nSteps = Length[Range[0, currentTime, dt]];
phaseData[[1 ;; nSteps]] =
Table[{i, normalize[x1[i]], x2[i]}, {i, 0, currentTime, dt}];
poincareMap =
makePoincreSection[q, torqueAmplitude, omega, phase, x10, x20,
isDamped, plotImageSize, axesLabel];
phasePlot =
makePhasePlot[phaseData, plotImageSize, axesLabel, nSteps];
bottomPlot = If[plotChoice == 1,
makeTimeSeriesPlot[phaseData, plotImageSize, nSteps,
currentTime],
makePowerSpectrumPlot[phaseData, plotImageSize, nSteps, omega,
spectrumScale, currentTime]
];
track = {Blue, PointSize[0.03],
Point[ { normalize@x1[currentTime], x2[currentTime]}]};
{nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3, maxSpeed}
];
(*---------------------------------------------------------------------------------*)
makeOneStep[
q_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),(*quality factor*)
torqueAmplitude_?(NumericQ[#] && Im[#] == 0 &&
Re[#] >= 0 &),(*driver amplitude*)
omega_?(NumericQ[#] &&
Im[#] == 0 &), (*driver frequency*)
phase_?(NumericQ[#] &&
Im[#] == 0 &), (*driver phase*)
x10_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for angle at time=0*)
x20_?(NumericQ[#] &&
Im[#] == 0 &), (*initial condition for speed at time=0*)
isDamped_?(Element[#, Booleans] &),
plotImageSize_List,
axesLabel_List,
plotChoice_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
pphaseData_,
spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
ccurrentTime_,
nnSteps_,
dt_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
pphasePlot_,
bbottomPlot_,
ttrack_,
maxRunTime_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
xx1_,
xx2_,
xx3_ ] :=
Module[{x1 = xx1, x2 = xx2, x3 = xx3, track = ttrack,
bottomPlot = bbottomPlot, phasePlot = pphasePlot,
nSteps = nnSteps, currentTime = ccurrentTime,
phaseData = pphaseData},
If[currentTime + dt <= maxRunTime,
(
nSteps += 1;
currentTime += dt;
{x1, x2, x3} =
solve[q, torqueAmplitude, omega, phase, currentTime,
currentTime + 1, x10, x20, isDamped, 7, 7];
phaseData[[nSteps]] = {currentTime, normalize@x1[currentTime],
x2[currentTime]};
phasePlot =
makePhasePlot[phaseData, plotImageSize, axesLabel, nSteps];
If[plotChoice == 1,
bottomPlot =
makeTimeSeriesPlot[phaseData, plotImageSize, nSteps,
currentTime],
bottomPlot =
makePowerSpectrumPlot[phaseData, plotImageSize, nSteps,
omega, spectrumScale, currentTime]
];
track = {Blue, PointSize[0.03],
Point[ { normalize@x1[currentTime], x2[currentTime] } ] };
)
];
{currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
x1, x2, x3}
];
(*---------------------------------------------------------------------------------*)
arcArrow[centre_, radius_, phi1_, phi2_, flip_, resolution_: 10] :=
Module[{i, data},
data = Table[
centre + radius {Cos[i], Sin[i]}, {i, phi1,
phi2, (phi2 - phi1)/Ceiling[(phi2 - phi1)/(Pi/resolution)]}
];
If[flip, data = Reverse[data]];
Arrow@data
]
)
]