Plot the results of a nonlinear system of equations

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;

• Post the Mathematica code (in code blocks) for your equations. See Markdown help – Bob Hanlon Jul 23 '17 at 16:54
• If you include the values of the parameters then readers would be able to try their solution and see that it works before posting it. – Bill Jul 23 '17 at 17:27
• Here are the values of the parameters that I'm sure of... alpha1=-24.85 alpha2=-66.319 alpha3=-9.2 alpha4=0 alpha5=-3.6 alpha6=-66.3 alpha7=0 the rest of the parameters should be set in order for the plot to take the desired form. @Bill – Soroush S. Jul 23 '17 at 18:42
• If somebody could just show me how to do it,I will change the parameters in order to reach the desired graph.Thank You very much – Soroush S. Jul 23 '17 at 18:51
• There are many ways to obtain the desired plot. I would use Solve to obtain expressions for the Sins of the two gammas in the first two equations, and then of the Coss off the same quantities from the last two equations. Square the respective results and add to eliminate the two gammas, leaving two polynomials in in alpha1 and alpha2. Then, use Solve again to obtain expressions for them, which can be plotted using Plot. – bbgodfrey Jul 23 '17 at 19:31

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.

• Thank you very much my friend.I'm grateful – Soroush S. Jul 24 '17 at 22:45

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