# Calculating the equilibria of a system of two coupled oscillators?

Equations: I have the following system of ODEs:

Phi1'[t]=omega1-a*Sin[Phi1[t]-Phi2[t]];
Phi2'[t]=-omega1-a*Sin[Phi2[t]-Phi1[t]];


Question: Is there a straightforward way of using Mathematica to calculate the equilibria of this system and determine their stability?

Comments: If this system did not entail the sine terms, I would simply set the equations equal to zero, solve for the equilibria, and determine the conditions under which the equilibria satisfy the Routh-Hurwitz criterion. But I am not sure how to deal with the sine terms. Any help would be greatly appreciated!

This has a conserved quantity since Phi1'[t]+Phi2'[t]==0. You can "algebraicize" by trig expanding and making everything a polynomial in sines and cosines of the phis. Then augment with trig identities (or use the Weierstrass substitution). I do the former below. Note that I handle the conservation also using a trig, and introducing a new constant.

exprs0 = {omega1 - a*Sin[phi1 - phi2], -omega1 - a*Sin[phi2 - phi1],
Cos[phi1 + phi2] - c};
exprs1 = TrigExpand[exprs0];
params = {omega1, a, c};
svars = Union[Cases[exprs1, Sin[a_], Infinity]];
cvars = Union[Cases[exprs1, Cos[a_], Infinity]];
augment = Map[#^2 + Cos[#[[1]]]^2 - 1 &, svars];
trigrule = {Cos -> cos, Sin -> sin};
exprs = Join[exprs1, augment] /. trigrule
tvars = Complement[Variables[exprs], params]

(* {omega1 - a cos[phi2] sin[phi1] + a cos[phi1] sin[phi2], -omega1 +
a cos[phi2] sin[phi1] - a cos[phi1] sin[phi2], -c +
cos[phi1] cos[phi2] - sin[phi1] sin[phi2], -1 + cos[phi1]^2 +
sin[phi1]^2, -1 + cos[phi2]^2 + sin[phi2]^2}

{cos[phi1], cos[phi2], sin[phi1], sin[phi2]} *)


I solve with a particular set of parameter values to find those trig values. The step of applying ArcSin to the sines is not necessarily correct since we are dealing with a multivalued inverse.

solns = SolveValues[(exprs == 0) /. Thread[params -> {1, 2, 1}],
tvars];
sines = solns[[All, 3 ;; 4]];
Map[ArcSin, sines, {2}] // N

(* Out[446]= {{-0.261799, 0.261799}, {0.261799, -0.261799}, {-1.309,
1.309}, {1.309, -1.309}} *)


Have a look at DSolveValue.

Here, I'm using initial conditions and some values for omega1 and a to plot a solution, but DSolveValue can indeed handle the fully-symbolic case too:

eqs = {
Phi1'[t] == omega1 - a*Sin[Phi1[t] - Phi2[t]],
Phi2'[t] == -omega1 - a*Sin[Phi2[t] - Phi1[t]]
};

ics = {Phi1[0] == 1, Phi2[0] == 0};

sols[t_] = DSolveValue[{eqs, ics}, {Phi1[t], Phi2[t]}, t]
Plot[Evaluate[sols[t] /. {omega1 -> 1, a -> 2}], {t, 0, 2 \[Pi]}]


• Thanks! I am going to give this a try. Dec 1, 2022 at 1:19
• Also, if by "equilibria" you mean the steady-state behavior - you can perhaps look at Limit or Asymptotic, e.g. Asymptotic[sols[t] /. {omega1 -> 1, a -> 2}, t -> \[Infinity]] Dec 1, 2022 at 1:22
• After further analysis, I found that the Limit function does work for this! Thank you! Do you by any chance know how one could determine whether the equilibrium/asymptote is stable? Would the Routh-Hurwitz criterion suffice for this? Dec 1, 2022 at 5:02