Plot the results of a nonlinear system of equations

Question

I have this system of four nonlinear equations:

the unknowns are a1, a2, gamma1 and gamma 2. all other parameters are known. What I want is to plot the a1 vs sigma2 as below:

I tried findinstance but it takes a lot of time to find the results for any value of sigma2, Let alone plotting a1 for different values of sigma2! Can anybody help me? Thanks in advance Here are the codes for the four equations and the corresponding values for the parameters:

8*Subscript[\[Omega], 1]*Subscript[\[Mu], 1]*Subscript[a, 1] - 
   Subscript[\[Alpha], 2]*Subscript[a, 1]^2*Subscript[a, 2]*
    Sin[Subscript[\[Gamma], 1]] - 
   4*Subscript[f, 1]*Sin[Subscript[\[Gamma], 2]] == 0; 
8*Subscript[\[Omega], 2]*Subscript[\[Mu], 2]*Subscript[a, 2] + 
   Subscript[\[Alpha], 5]*Subscript[a, 1]^3*
    Sin[Subscript[\[Gamma], 1]] == 0; 
8*Subscript[\[Omega], 1]*Subscript[a, 1]*
    Subscript[\[Sigma], 
     2] + (3*Subscript[\[Alpha], 1]*Subscript[a, 1]^2 + 
      2*Subscript[\[Alpha], 3]*Subscript[a, 2]^2)*Subscript[a, 1] + 
       Subscript[\[Alpha], 2]*Subscript[a, 1]^2*Subscript[a, 2]*
    Cos[Subscript[\[Gamma], 1]] + 
   4*Subscript[f, 1]*Cos[Subscript[\[Gamma], 2]] == 0; 
8*Subscript[\[Omega], 2]*
    Subscript[a, 
     2]*(3*Subscript[\[Sigma], 2] - 
      Subscript[\[Sigma], 1]) + (3*Subscript[\[Alpha], 8]*
       Subscript[a, 2]^2 + 
      2*Subscript[\[Alpha], 6]*Subscript[a, 1]^2)*Subscript[a, 2] + 
       Subscript[\[Alpha], 5]*Subscript[a, 1]^3*
    Cos[Subscript[\[Gamma], 1]] == 0; 

Also the values for the parameters are:

Subscript[\[Omega], 1] = 10; 
Subscript[\[Mu], 1] = 1; 
Subscript[\[Alpha], 2] = -66.319; 
Subscript[f, 1] = 1; 
Subscript[\[Omega], 2] = 30; 
Subscript[\[Mu], 2] = 1; 
Subscript[\[Alpha], 5] = -3.6; 
Subscript[\[Alpha], 1] = -24.85; 
Subscript[\[Alpha], 3] = -9.2; 
Subscript[\[Alpha], 6] = -66.3; 
Subscript[\[Alpha], 8] = -345; 
Subscript[\[Sigma], 1] = 0.08; 
f1 = 0; 

Answer 1

Generally, avoiding Subscript is a good idea. Therefore, with the four equations in the question designated as eq, transform them to

eq = eq /. Subscript[z_, n_] -> z[n]
(* {-4 f[1] Sin[γ[2]] - a[1]^2 a[2] Sin[γ[1]] α[2] + 8 a[1] μ[1] ω[1] == 0, 
    a[1]^3 Sin[γ[1]] α[5] + 8 a[2] μ[2] ω[2] == 0, 
    4 Cos[γ[2]] f[1] + a[1]^2 a[2] Cos[γ[1]] α[2] + a[1] (3 a[1]^2 α[1] + 2 a[2]^2 α[3]) + 
        8 a[1] σ[2] ω[1] == 0, 
    a[1]^3 Cos[γ[1]] α[5] + a[2] (2 a[1]^2 α[6] + 3 a[2]^2 α[8]) + 
        8 a[2] (-σ[1] + 3 σ[2]) ω[2] == 0} *)

Similarly, the parameter values, rationalized and recast as rules for convenience, are designate rules and given by

{ω[1] -> 10, μ[1] -> 1, α[2] -> -(66319/1000), f[1] -> 1, ω[2] -> 30, μ[2] -> 1, 
 α[5] -> -(18/5), α[1] -> -(497/20), α[3] -> -(46/5), α[6] -> -(663/10), 
 α[8] -> -345, σ[1] -> 2/25}

Then, γ[1] and γ[2] are eliminated by

ss = Flatten@Solve[eq[[1 ;; 2]] /. rules, {Sin[γ[1]], Sin[γ[2]]}] // Factor;
sc = Flatten@Solve[eq[[3 ;; 4]] /. rules, {Cos[γ[1]], Cos[γ[2]]}] // Factor;
eqs = {Sin[γ[1]]^2 + Cos[γ[1]]^2 == 1, Sin[γ[2]]^2 + Cos[γ[2]]^2 == 1} /. ss /. sc
(* {(40000 a[2]^2)/(9 a[1]^6) + (a[2]^2 (32 + 221 a[1]^2 + 1725 a[2]^2 - 1200 σ[2])^2)
        /(36 a[1]^6) == 1,
    (1200 a[1]^2 + 66319 a[2]^2)^2/(3600 a[1]^2) + (447300 a[1]^4 - 2122208 a[2]^2 - 
         14546099 a[1]^2 a[2]^2 - 114400275 a[2]^4 - 480000 a[1]^2 σ[2] + 
         79582800 a[2]^2 σ[2])^2/(576000000 a[1]^2) == 1} *)

Because these equations represent a 32-order polynomial system, obtaining numerical rather than symbolic solutions is the more practical option. A scan of the range of {σ[2], -10, 35} shown in the figure in the question yields,

Table[tem = NSolve[eqs, {a[1], a[2]}, Reals, WorkingPrecision -> 30]; 
    If[tem == {}, Nothing, {σ[2], #} & /@ Union[a[1] /. tem, 
    SameTest -> (Abs[#1 - #2] < 10^-5 &)]], {σ[2], -10, 35}]
(* {{{-1, -0.0353347773347261876634968709486}, {-1, .0353347773347261876634968709486}}, 
    {{0, 0.049999860456465014315945496498}, {0, .049999860456465014315945496498}}, 
    {{1, -0.0353759608868205937784722694063}, {1, 0.0353759608868205937784722694063}}} *)

Evidently, real solutions occur only near σ[2] == 0. Providing greater resolution there (i.e., {σ[2], -3/2, 3/2, 1/20}) yields,

ListPlot[%, AxesLabel -> {σ[2], a[1]}, PlotStyle -> Blue, 
    ImageSize -> Large, LabelStyle -> Directive[12, Black, Bold]]

That this result differs markedly from those in the figure in the question should not be surprising. The value of f[1] specified in the question differs from those in the plot by two orders of magnitude or more. Increasing f[1] to 100 yields a curve qualitatively similar to some of those in the question. Perhaps, some of the other parameter values used to obtain that plot also differ from those specified in the question. In any case, this seems like a credible result.

Answer 2

Here's a numerical approach using pseudo-arclength continuation using a function I wrote here and borrowing @bbgodgrey's de-subscripted equations.

First, define TrackRootPAL from here. Then set up the equations and parameter values:

eq = {-4 f[1] Sin[γ[2]] - a[1]^2 a[2] Sin[γ[1]] α[2] + 8 a[1] μ[1] ω[1],
  a[1]^3 Sin[γ[1]] α[5] + 8 a[2] μ[2] ω[2],
  4 Cos[γ[2]] f[1] + a[1]^2 a[2] Cos[γ[1]] α[2] + a[1] (3 a[1]^2 α[1] + 2 a[2]^2 α[3]) + 8 a[1] σ[2] ω[1], 
  a[1]^3 Cos[γ[1]] α[5] + a[2] (2 a[1]^2 α[6] + 3 a[2]^2 α[8]) + 8 a[2] (-σ[1] + 3 σ[2]) ω[2]};

rules = {ω[1] -> 10, μ[1] -> 1, α[2] -> -(66319/1000), f[1] -> 100, ω[2] -> 30, μ[2] -> 1, α[5] -> -(18/5), α[1] -> -(497/20), α[3] -> -(46/5), α[6] -> -(663/10), α[8] -> -345, σ[1] -> 2/25};

Then track the solution using TrackRootPAL[eqns, unks, {par, parmin, parmax}, initguesspar, initguess]:

tr = TrackRootPAL[eq /. rules, {a[1], a[2], γ[1], γ[2]}, {σ[2], -10, 35},
  0, {1.5, 0, 0, 0}];

(* 3 solution branches returned as InterpolatingFunctions *)

Plot[Evaluate[a[1][σ[2]] /. tr], {σ[2], -10, 35}]