# Trying to solve a system of non-linear equations

I am trying to use Mathematica to find a solution for the following system of equations:

q=a0*(w/p*e)^a1
e=b0*q^b1*(w/p)^(-b2)
p=c0*w^c1*q^c2
w=d0*p^d1*e^d2


where a0, a1, b0, b1 etc. are parameters and q, e, p and w are endogenous variables. I would like to end up with expressions for q, e, p and w which are only in terms of the parameters. Trying Solve on the above only returns expressions for p and q in terms of parameters, w and e, along with the warnings:

Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information

and

Equations may not give solutions for all "solve" variables

Trying Reduce instead results in Mathematica computing for an extended period until I eventually abort.

The commands I use are as follows:

eqn = {
q == Subscript[a, 0]*(w/p*e)^Subscript[a, 1]
, e == Subscript[b, 0]*q^Subscript[b, 1]*(w/p)^(-Subscript[b, 2])
, p == Subscript[c, 0]*w^Subscript[c, 1]*q^Subscript[c, 2]
, w == Subscript[d, 0]*p^Subscript[d, 1]*e^Subscript[d, 2]
}

Solve[eqn, {w,e,p,q}]


and

Reduce[eqn, {w,e,p,q}]


I have also tried to clarify that all the parameters should be >0 by adding the qualifications &&Subscript[a, 0]>0 && Subscript[a, 1]>0 etc. to the definition of eqn but this results in Mathematica seemingly being unable to do anything with the expressions in eqn.

I am not sure whether this is just a general problem with my equations or whether I am doing something wrong since I am not very experienced in Mathematica. Any help would be much appreciated.

• In the definition of eqn you are using emp rather than e as shown in the first set. In the call to Solve and Reduce you are using e as a variable. You need to either use emp or e in both places, the definition of as well as in Solve. – Jack LaVigne Aug 27 '17 at 14:00
• Thanks, that is an error which slipped in when posting the question - it's not present in the code that I'm running so can't be the source of the problem. – SDR91 Aug 27 '17 at 14:45

Not sure if this is the best way to go about it but you might take logs of both sides, using the assumptions to simplify so that it becomes a linear system.

eqns = {q == a0*(w/p*e)^a1,
e == b0*q^b1*(w/p)^(-b2),
p == c0*w^c1*q^c2,
w == d0*p^d1*e^d2};
logeqns = Map[Log, eqns, {2}]

(* Out= {Log[q] == Log[a0 ((e w)/p)^a1],
Log[e] == Log[b0 q^b1 (w/p)^-b2], Log[p] == Log[c0 q^c2 w^c1],
Log[w] == Log[d0 e^d2 p^d1]} *)

logeqns2 =
PowerExpand[logeqns,
Assumptions -> {a0 > 0, a1 > 0, b0 > 0, b1 > 0, b2 > 0, c0 > 0,
c1 > 0, c2 > 0, d0 > 0, d1 > 0, d2 > 0, q > 0, e > 0, p > 0,
w > 0}]

(* Out= {Log[q] == Log[a0] + a1 (Log[e] - Log[p] + Log[w]),
Log[e] == Log[b0] + b1 Log[q] - b2 (-Log[p] + Log[w]),
Log[p] == Log[c0] + c2 Log[q] + c1 Log[w],
Log[w] == Log[d0] + d2 Log[e] + d1 Log[p]} *)


Make new variables for the logs and solve it.

logeqns2 =
PowerExpand[logeqns,
Assumptions -> {a0 > 0, a1 > 0, b0 > 0, b1 > 0, b2 > 0, c0 > 0,
c1 > 0, c2 > 0, d0 > 0, d1 > 0, d2 > 0, q > 0, e > 0, p > 0,
w > 0}]

(* Out= {Log[q] == Log[a0] + a1 (Log[e] - Log[p] + Log[w]),
Log[e] == Log[b0] + b1 Log[q] - b2 (-Log[p] + Log[w]),
Log[p] == Log[c0] + c2 Log[q] + c1 Log[w],
Log[w] == Log[d0] + d2 Log[e] + d1 Log[p]} *)

lgeqns = logeqns2 /. Log[a_] :> lg[a];
lgsolns = Solve[lgeqns, {lg[q], lg[e], lg[p], lg[w]}];
solns = Exp[{lg[q], lg[e], lg[p], lg[w]} /. lgSolns /. lg -> Log]

(* Out= {{E^(-((-Log[a0] + c1 d1 Log[a0] - b2 d2 Log[a0] +
b2 c1 d2 Log[a0] - a1 Log[b0] + a1 c1 d1 Log[b0] - a1 d2 Log[b0] +
a1 c1 d2 Log[b0] + a1 Log[c0] - a1 b2 Log[c0] - a1 d1 Log[c0] +
a1 b2 d1 Log[c0] - a1 Log[d0] + a1 b2 Log[d0] + a1 c1 Log[d0] -
a1 b2 c1 Log[d0])/(
1 - a1 b1 + a1 c2 - a1 b2 c2 - c1 d1 + a1 b1 c1 d1 - a1 c2 d1 +
a1 b2 c2 d1 - a1 b1 d2 + b2 d2 + a1 b1 c1 d2 - b2 c1 d2))),
E^(-((-b1 Log[a0] - b2 c2 Log[a0] + b1 c1 d1 Log[a0] +
b2 c2 d1 Log[a0] - Log[b0] - a1 c2 Log[b0] + c1 d1 Log[b0] +
a1 c2 d1 Log[b0] + a1 b1 Log[c0] - b2 Log[c0] - a1 b1 d1 Log[c0] +
b2 d1 Log[c0] - a1 b1 Log[d0] + b2 Log[d0] + a1 b1 c1 Log[d0] -
b2 c1 Log[d0])/(
1 - a1 b1 + a1 c2 - a1 b2 c2 - c1 d1 + a1 b1 c1 d1 - a1 c2 d1 +
a1 b2 c2 d1 - a1 b1 d2 + b2 d2 + a1 b1 c1 d2 - b2 c1 d2))),
E^(-((-c2 Log[a0] - b1 c1 d2 Log[a0] - b2 c2 d2 Log[a0] -
a1 c2 Log[b0] - c1 d2 Log[b0] - a1 c2 d2 Log[b0] - Log[c0] +
a1 b1 Log[c0] + a1 b1 d2 Log[c0] - b2 d2 Log[c0] - c1 Log[d0] +
a1 b1 c1 Log[d0] - a1 c2 Log[d0] + a1 b2 c2 Log[d0])/(
1 - a1 b1 + a1 c2 - a1 b2 c2 - c1 d1 + a1 b1 c1 d1 - a1 c2 d1 +
a1 b2 c2 d1 - a1 b1 d2 + b2 d2 + a1 b1 c1 d2 - b2 c1 d2))),
E^(-((-c2 d1 Log[a0] - b1 d2 Log[a0] - b2 c2 d2 Log[a0] -
a1 c2 d1 Log[b0] - d2 Log[b0] - a1 c2 d2 Log[b0] - d1 Log[c0] +
a1 b1 d1 Log[c0] + a1 b1 d2 Log[c0] - b2 d2 Log[c0] - Log[d0] +
a1 b1 Log[d0] - a1 c2 Log[d0] + a1 b2 c2 Log[d0])/(
1 - a1 b1 + a1 c2 - a1 b2 c2 - c1 d1 + a1 b1 c1 d1 - a1 c2 d1 +
a1 b2 c2 d1 - a1 b1 d2 + b2 d2 + a1 b1 c1 d2 - b2 c1 d2)))}} *)


Not very pretty, due to that exponential denominator. Maybe could be simplified some, I'm not sure.

• Thanks, this does work but as you say the solutions are not pretty and aren't very useful for further processing. Perhaps I'll try to rewrite the system to get rid of at least some of the nonlinearities. – SDR91 Aug 27 '17 at 15:31