Solving Non-linear System of Equations

Question

I need to get a symbolic solution to a system of non-linear equations (equation 8 is non-linear), but Mathematica is running for ever and gives no answer. Is there another way to get the solution? My codes are:

eq1 = g1 - g2 - g3 == 0
eq2 = (a1 + a2) (g1 - gA) - a2 g2 == 0
eq3 = g2 - ((b1 + b2)/b1) gB - ((1. - b1 - b2)/b1) g3 == 0
eq4 = g3 - r - b2 g2 (p - 1.) == 0
eq5 = r - ρ - g1 == 0
eq6 = gA - ηA LA == 0
eq7 = gB - ηB (L - Ly - Lq - LA) - f g2 == 0
eq8 = (p - 1.) b1 Ly/ a1 Lq - p/ a2 == 0
eq9 = a1 (a1 + a2) r - a1 a2 g2 - ηA (a1 + a2)^2. (1 - a1 - a2) Ly == 0
eq10 = b1 (b1 + b2) r - b1 g3 - b2 b1 g2 - ηB (b1 + b2)^2. (1. - b1 - b2) Lq == 0
aa = Solve[{eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10}, {Ly, Lq, LA, gA, gB, r,g1,g2,g3,p}]

Thanks, Hossein

Answer 1

These equations are solvable, but the process is exceptionally slow and the output huge. To help Solve, replace the approximate real numbers 1. and 2. by 1 and 2. Then, solve the nine linear equations in terms of p and the parameters.

nine = Simplify[Solve[{eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq9, eq10}, {Ly, Lq, LA, 
      gA, gB, r, g1, g2, g3}][[1]]];

and substitute the result into eq8.

e8 = Simplify[eq8 /. nine]
(* (b1^3*(b1 + b2)*(-1 + p)*(-1 + b2 - b2*p)*((-1 + a2)*a2^2*((-1 + b1 + b2)*L*(1 + b2*(-1 + p))*ηA*ηB + 
      ((-1 + b1 + b2)*(1 - f + b2*(-1 + p))*ηA + b2*(1 + b1*(-1 + p) + b2*(-1 + p))*ηB)*ρ) + 
    a1^3*((-1 + b1 + b2)*L*(1 + b2*(-1 + p))*ηA*ηB + ((-1 + b1 + b2)*(1 - f + b2*(-1 + p))*ηA + 
        (-1 + 2*b2 + b1*(1 + b2*(-1 + p)) + b2^2*(-1 + p))*ηB)*ρ) + a1^2*((-1 + 3*a2)*(-1 + b1 + b2)*L*(1 + b2*(-1 + p))*ηA*ηB + 
      ((-1 + 3*a2)*(-1 + b1 + b2)*(1 - f + b2*(-1 + p))*ηA + (1 - b1 - 2*b2 + b1*b2 + b2^2 + 
          a2*(-2 + 5*b2 + b1*(2 + 3*b2*(-1 + p)) + 3*b2^2*(-1 + p)) - b2*p - b1*b2*p - b2^2*p)*ηB)*ρ) + 
    a1*a2*((-2 + 3*a2)*(-1 + b1 + b2)*L*(1 + b2*(-1 + p))*ηA*ηB + ((-2 + 3*a2)*(-1 + b1 + b2)*(1 - f + b2*(-1 + p))*ηA + 
        (a2*(-1 + 4*b2 + b1*(1 + 3*b2*(-1 + p)) + 3*b2^2*(-1 + p)) - b2*(2 + 2*b1*(-1 + p) + 2*b2*(-1 + p) + p))*ηB)*ρ))*
   (-((a2*(-1 + b1 + b2)*(b1^2*(ηA + ηB) + 2*b1*((-1 + b2)*ηA + b2*ηB) + b2*((-1 + b2)*ηA + b2*ηB))*η)/b1) + 
    ((a1 + a2)*ηA - (a1 + a2)*b2*ηA + ((a1 + a2)*(-1 + b1 + b2)^2*(b1 + b2)*ηA)/b1 + 
      ((-1 + b1 + b2)*(b1 + b2)*(a1*(b1*(-1 + f) + b2*f)*ηA + a2*(b1*((-1 + f)*ηA + ηB) + b2*(f*ηA + ηB))))/b1)*ρ + 
    (-1 + b2 - b2*p)*(-((-a1 - a2)*(1 + ((-1 + b1 + b2)^2*(b1 + b2))/b1)*ηA*ρ) + ((a1 + a2)*(-1 + b1 + b2)*(b1 + b2)^2*ηB*(L*ηA + ρ))/b1)))/
  ((a1 + a2)^2*(-1 + b1 + b2)*(1 + b2*(-1 + p))^3*ηA*ηB*(-(a1*(-1 + b1 + b2)*(b1 + b2)^3*ηB) + 
    (-1 + a1 + a2)*(a1 + a2)*(b1^2 + (-1 + b2)*b2 + b1*(-1 + 2*b2))*(b1^2*(ηA + ηB) + 2*b1*((-1 + b2)*ηA + b2*ηB) + 
      b2*((-1 + b2)*ηA + b2*ηB)))*(a1^2*(b1^2*(ηA + ηB) + 2*b1*((-1 + b2)*ηA + b2*ηB) + b2*((-1 + b2)*ηA + b2*ηB)) + 
    (-1 + a2)*a2*(b1^2*(ηA + ηB) + 2*b1*((-1 + b2)*ηA + b2*ηB) + b2*((-1 + b2)*ηA + b2*ηB)) + 
    a1*(b1^2*((-1 + 2*a2)*ηA + 2*(-1 + a2)*ηB) + 2*b1*((-1 + 2*a2)*(-1 + b2)*ηA + 2*(-1 + a2)*b2*ηB) + 
      b2*((-1 + 2*a2)*(-1 + b2)*ηA + 2*(-1 + a2)*b2*ηB)))) == p/a2 *)

Finally, obtain the three roots of e8

Solve[e8, p]

any one of which Mathematica describes as "large output" and does not display. I imagine that the computation would be much faster and the output much more compact, if the parameters were replaced by numbers.

Alternative Approach

Because the nine linear equations are sparse, as can be seen from

Normal[CoefficientArrays[{eq1, eq4, eq5}, {g1, r, g3}][[2]]]

it is possible to eliminate all variables but p more or less one at a time. {eq1, eq4, eq5} yield

{g1 -> r - ρ, g2 -> -(ρ/(1 - b2 + b2 p)), g3 -> r - ρ + ρ/(1 - b2 + b2 p)}

Subsequent back substitutions and eliminations then yield, in turn,

gA -> LA ηA
gB -> L ηB - (LA ηB + Lq ηB + Ly ηB + (f ρ)/(1 + b2 (-1 + p)))
r -> -((-1 + a1 + a2) (a1 + a2)^2 Ly ηA + (a1 a2 ρ)/(1 + b2 (-1 + p)))/(a1 (a1 + a2))
LA -> (a1 Ly ηA - a1^2 Ly ηA + a2 Ly ηA - 
  2 a1 a2 Ly ηA - a2^2 Ly ηA - a1 ρ)/(a1 ηA)
Lq -> (-((b1 b2 ρ)/(1 + b2 (-1 + p))) - (b1 (b1 + b2) (-(-1 + a1 + a2) (a1 + a2)^2 Ly ηA
  - (a1 a2 ρ)/(1 + b2 (-1 + p))))/(a1 (a1 + a2)) + b1 (-ρ + ρ/
  (1 + b2 (-1 + p)) - ((-1 + a1 + a2) (a1 + a2)^2 Ly ηA + (a1 a2 ρ)/
  (1 + b2 (-1 + p)))/(a1 (a1 + a2))))/((-1 + b1 + b2) (b1 + b2)^2 ηB)
Ly -> (-(((b1 + b2) L ηB)/b1) + ((1 - b1 - b2) ρ)/b1 - ρ/((-1 + b1 + b2) (b1 + b2)) - 
  ρ/(1 + b2 (-1 + p)) - ((1 - b1 - b2) ρ)/(b1 (1 + b2 (-1 + p))) + (a2 (1 - b1 - b2) ρ)/
  ((a1 + a2) b1 (1 + b2 (-1 + p))) + (a2 ρ)/((a1 + a2) (-1 + b1 + b2) (1 + b2 (-1 + p))) + 
  ρ/((-1 + b1 + b2) (b1 + b2) (1 + b2 (-1 + p))) - (a2 ρ)/((a1 + a2) (-1 + b1 + b2) 
  (b1 + b2) (1 + b2 (-1 + p))) - (b2 ρ)/((-1 + b1 + b2) (b1 + b2) (1 + b2 (-1 + p))) + 
  ((b1 + b2) f ρ)/(b1 (1 + b2 (-1 + p))) - ((b1 + b2) ηB ρ)/(b1 ηA))/
  (-(((-1 + a1 + a2) (a1 + a2) (1 - b1 - b2) ηA)/(a1 b1)) - ((-1 + a1 + a2) (a1 + a2) ηA)/
  (a1 (-1 + b1 + b2)) + ((-1 + a1 + a2) (a1 + a2) ηA)/(a1 (-1 + b1 + b2) (b1 + b2)) -
  ((b1 + b2) ηB)/b1 + (a1 (b1 + b2) ηB)/b1 + ((-1 + a2) a2 (b1 + b2) ηB)/
  (a1 b1) + ((-1 + 2 a2) (b1 + b2) ηB)/b1)

This leaves eq8, which is too long to be reproduced here, which can be solved as in the earlier approach.

Answer 2

We can rationalize to make everything exact, clear denominators, and compute a Groebner basis (slow, but it works) over the rational function field involving the parameters (non-variables, that is).

polys = Numerator[
   Together[
    Rationalize[
     Map[First, {eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, 
       eq10}]]]];
vars = {Ly, Lq, LA, gA, gB, r, g1, g2, g3, p};

Timing[
 gb = GroebnerBasis[polys, vars, 
    CoefficientDomain -> RationalFunctions];]

(* Out[190]= {726.063, Null} *)

It's not small...

In[191]:= LeafCount[gb]

(* Out[191]= 4567713 *)

Structure is convenient however. First polynomial is cubic in p, the rest are each linear in a unique variable and quaratic in p. So each of the three solutions for p will produce one value for each remaining variable.

In[192]:= Length[gb]

(* Out[192]= 10 *)

In[198]:= Map[Intersection[Variables[#], vars] &, gb]

(* Out[198]= {{p}, {g3, p}, {g2, p}, {g1, p}, {p, r}, {gB, p}, {gA, 
  p}, {LA, p}, {Lq, p}, {Ly, p}} *)

In[199]:= Exponent[gb[[1]], p]

(* Out[199]= 3 *)

All this aside, I would actually recommend a numeric approach wherein NSolve is invoked for inputs where the parameters all have numeric instantiations. This is less prone to numeric issues resulting from trying to numerically instantiate within a massive symbolic solution of the sort that Solve might provide from the Groebner basis.

--- edit ---

I should mention that adding Sort->True to the GroebnerBasis computation brings it down to 160 seconds on the same machine, with a result that is around 1/9 the size (measured by LeafCount). So it's not entirely hopeless to do this symbolically.

--- end edit ---