Trying to solve a system of non-linear equations

Question

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.

Answer 1

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[231]= {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[235]= {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[235]= {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[252]= {{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.