0
$\begingroup$

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!

$\endgroup$

2 Answers 2

1
$\begingroup$

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}} *)
$\endgroup$
2
$\begingroup$

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]}]

enter image description here

$\endgroup$
3
  • 1
    $\begingroup$ Thanks! I am going to give this a try. $\endgroup$
    – MathIsHard
    Dec 1, 2022 at 1:19
  • 1
    $\begingroup$ 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]] $\endgroup$ Dec 1, 2022 at 1:22
  • $\begingroup$ 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? $\endgroup$
    – MathIsHard
    Dec 1, 2022 at 5:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.