# Solving Non-linear System of Equations

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

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Jan 20, 2015 at 0:39
• You have more variables than equations. Jan 20, 2015 at 1:10
• There are 10 variables and 10 equations. Jan 20, 2015 at 1:32
• These are my parameters: {a1, a2, b1, b2,[Rho], [Eta]A, [Eta]B, f, L} and these are my 10 variables: {Ly, Lq, LA, gA, gB, r, g1, g2, g3, p}, I need the solution: variables as functions of the parameters. Jan 20, 2015 at 2:13
• This subsystem seems to be inconsistent: Solve[{eq1, eq4, eq5}, {r, g1, g3}] Jan 20, 2015 at 3:24

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.

• I think the preprocessing approach of dealing with the linear equations first is probably best. Definitely better than the brute-force method I showed. Aug 28, 2016 at 15:28

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