I'm looking at the behavior of a system of ODEs as a function of a parameter.
This is a system of phase oscillators, with the dynamics (simplified for clarity of reading, proper code farther below):
Phi'[m][t] == Freq + Sin[Phi[m + 1][t] - Phi[m][t] - p]
Thanks to helpful answers here (Plotting interpolating functions from ParametricNDSolveValue?), I understand how to plot the difference between two dynamic variables (for example, dPhi = Phi[m+1]-Phi[m]
) for different values of the parameter, p.
I'm still struggling to plot dPhi as a function of parameter, p. This system has a transient before reaching a steady state, and I'd like to plot this steady state phase difference versus the parameter p, either as a Plot or a ListPlot.
(*system parameters*)
a = RotateLeft[IdentityMatrix[4], {1, 0}];
NatFreq = {1, .98, 1.02, 1.01};
InitCond[n_] := Array[Pi (# - 1)/(4) &, n]
(*ODEs and initial conditions*)
deqns =
Table[
Phi[m]'[t] ==
NatFreq[[m]] + Sum[a[[m,j]]*Sin[Phi[j][t] - Phi[m][t] - p],{j,4}],
{m, 1, 4}]
ics = Table[Phi[m][0] == InitCond[4][[m]], {m, 1, 4}];
vars = Table[Phi[m], {m, 4}];
pVals = {0, Pi/2, Pi, 3 Pi/2};
(*compute the solution*)
pfun = ParametricNDSolveValue[{deqns, ics}, vars, {t, 0, 100}, {p}];
pPhi[n : (1 | 2 | 3 | 4), p_] := pfun[p][[n]]
With this, I can make a lovely plot of the transient of the phase differences for several values of p:
d41 = Mod[
Evaluate[
Table[pPhi[4,p][t], {p, pVals}] -
Table[pPhi[1,p][t], {p, pVals}]
],
2*Pi];
Plot[d41, {t, 0, 100},
PlotRange -> {All, {0, 2*Pi}},
PlotStyle -> {Black, Red, Green, Blue},
AxesLabel -> {time, Phi4 - Phi1},
PlotLegends -> pVals]
What I haven't been able to figure out is how to plot the final values of these transients against the parameter p rather than against time. This was my best guess:
T=50;
Plot[Evaluate[pPhi[2,p][t] /. t -> T-pPhi[1,p][t] /. t -> T], {p, 0, 2*Pi},
PlotRange -> {{0, 2*Pi}, {0, 2*Pi}]
Everything I try produces many errors. This seems like it may be very simple, but I haven't been able to work my way through it. (Please also note that it must be NDSolve related and not DSolve; I use Sine in this is example, but I sometimes use a longer Fourier series.) Thanks very much.
Plot[Mod[pPhi[4, p][100] - pPhi[1, p][100], 2 \[Pi]], {p, 0, 2 \[Pi]}]
. $\endgroup$