1

I am trying to solve a system of four equations in four variables. I have read a number of threads on similar issues and tried to follow the suggestions. But I think it is a bit messy because of the logs and cross products here. This is the exact system:

7*w = (7*w+5*x+2*y+z) * ( 0.76 + 0.12*Log[w] -0.08*Log[x] -0.03*Log[y] -0.07*Log[7*w+5*x + 2*y + z]),
5*x = (7*w+5*x+2*y+z) * ( 0.84 - 0.08*Log[w] +0.11*Log[x] -0.02*Log[y] -0.08*Log[7*w+5*x + 2*y + z]),
2*y = (7*w+5*x+2*y+z) * (-0.45 - 0.03*Log[w] -0.02*Log[x] +0.05*Log[y] +0.12*Log[7*w+5*x + 2*y + z]),
1*z = (7*w+5*x+2*y+z) * (-0.16 + 0*Log[w] - 0*Log[x] - 0*Log[y] + 0.03*Log[7*w+5*x + 2*y + z])

(FYI-I am a young economist, and this is an extension of a consumer demand system.) Theoretically, we know that there exists a unique solution to this system that is positive.

Trys

  • Solve & NSolve : As there should be a solution I tried these but neither works. I guess that the system has too many logs to handle.

  • FindRoot : I started with an initial value of (14,15,10,100) which I get from my data. FindRoot returns the last value (which does not satisfy my system) and the following message.

    FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable.....

I tried different initial values, including the value returned by FindRoot. I tried to analyze the pattern of the solution value at each step. I didn’t see any pattern but noticed that the z values become negative early in the process. So I put bounds on the values. This just stops the code at the minimum value of 0.1. I also tried an exponential system instead of log, same issues.

 Reap[FindRoot[{
 7*w==(7*w+5*x + 2*y + z)*(0.76 + 0.12*Log[w] -0.08*Log[x] -0.03*Log[y] -0.07*Log[7*w+5*x + 2*y + z]),
 5*x==(7*w+5*x + 2*y + z)*(0.84  -0.08*Log[w] +0.11*Log[x] -0.02*Log[y] -0.08*Log[7*w+5*x + 2*y + z]),
 2*y==(7*w+5*x + 2*y + z)*(-0.45 - 0.03*Log[w] -0.02*Log[x] +0.05*Log[y] +0.12*Log[7*w+5*x + 2*y + z]),
 z==(7*w+5*x + 2*y + z)*(-0.16 + 0*Log[w] -0*Log[x] -0*Log[y] +0.03*Log[7*w+5*x + 2*y + z])},
      {{w,14,0.1,500},{x,15,0.1,500},{y,10,0.1,500},        
       {z,100,0.1,500}},EvaluationMonitor:>Sow[{w,x,y,z}] ]]
  • FindMinimum : As we can write this problem as a minimization problem, I tried this (following the suggestion here). The value returned did not converge the system or equations to zero. I tried with only the first two equations, and that sort of converged to zero.

Hope this is engaging enough for the experts here! Any ideas how I should find the solution or why can’t I? It’s the first time I am using Mathematica, and unfortunately the first time I am empirically solving a system/optimizing! Thanks a lot.

{g1,g2,g3, g4}={7*w - (7*w+5*x+2*y+z)* (0.76+0.12*Log[w]-0.08*Log[x]-0.03*Log[y] -0.07*Log[7*w+5*x+2*y+z]),5*x - (7*w+5*x+2*y+z)*(0.84-0.08*Log[w]+0.11*Log[x]-0.02*Log[y] -0.08*Log[7*w+5*x+2*y+z]),2*y - (7*w+5*x+2*y+z)*(-0.45-0.03*Log[w]-0.02*Log[x]+0.05*Log[y]+0.12*Log[7*w+5*x+2*y+z]), 1*z - (7*w+5*x+2*y+z)*(-0.16+0*Log[w]-0*Log[x]-0*Log[y]+0.03*Log[7*w+5*x+2*y+z])};subdomain=0<w<100 &&0<x<100 && 0<y<100 && 0<z<100;res=FindMinimum[{Total[{g1,g2,g3, g4}^2],subdomain},{w,x,y,z},AccuracyGoal->5]{g1,g2,g3,g4}/.res[[2]]
agentp
  • 6,849
  • 2
  • 19
  • 37
Divergent-Economist
  • 13
  • 3

2 Answers2

0

I don't have access to Mathematica, I put your equations into AMPL which is free for students. Here is what I did:

var w := 14 >= 0;
var x := 15 >= 0;
var y := 10 >= 0;
var z := 100 >= 0;

eq1: 7*w = (7*w+5*x+2*y+z) * ( 0.76 + 0.12*log(w) -0.08*log(x) -0.03*log(y) -0.07*log(7*w+5*x + 2*y + z));
eq2: 5*x = (7*w+5*x+2*y+z) * ( 0.84 - 0.08*log(w) +0.11*log(x) -0.02*log(y) -0.08*log(7*w+5*x + 2*y + z));
eq3: 2*y = (7*w+5*x+2*y+z) * (-0.45 - 0.03*log(w) -0.02*log(x) +0.05*log(y) +0.12*log(7*w+5*x + 2*y + z));
eq4: 1*z = (7*w+5*x+2*y+z) * (-0.16 +    0*log(w) -   0*log(x) -   0*log(y) +0.03*log(7*w+5*x + 2*y + z));

option show_stats 1;
option presolve 10;
option solver "/home/ali/ampl/ipopt"; # put your path here

option seed 1731;

# Initial solve
solve;
display w, x, y, z;

# Multistart
for {1..10} {
    for {j in 1.._snvars}
        let _svar[j] := Uniform(1, 50);
    solve;
    if (solve_result_num < 200) then {
      display w, x, y, z;
    }
}

If I only require that the variables are nonnegative, I get rubbish, for example

w = 2.39266e-11
x = 6.62678e-11
y = 1.57043e-24
z = 7.0842e-10

or 

w = 1.09972e-12
x = 9.77807e-11
y = 3.36229e-21
z = 1.85441e-09

Numerically, these are indeed solutions, they satisfy the equations to a fairly high precision, although I am pretty sure it's not what you are looking for. This indicates issues with your model.

If I increase the lower bounds of the variables a bit:

var w := 14 >= 0.1;
var x := 15 >= 0.1;
var y := 10 >= 0.1;
var z := 100 >= 0.01;

I get, even with multistart, Ipopt 3.11.6: Converged to a locally infeasible point. Problem may be infeasible. This again indicates issues with your model equations.

I am afraid you will have to revise your model.

This won't fix the issues with your model equations but I would introduce new variables: a=log(w), b=log(x), c=log(y), d=log(z). Then w=exp(a) and so on. It has the advantage that the function evaluations won't fail due to negative arguments for the logarithm function.

I would probably also introduce a new variable for (7*w+5*x+2*y+z) just to make the equations more compact.

Neither of these new variables will solve the above issues with your model equations.

If it is really your first time using Mathematica, you might be better off with AMPL and IPOPT; these tools are custom-tailored for solving equations and optimization problems. I suggest you use the AMPL mailing list if you have question and not Stackoverflow; you will get better answers on the mailing list.

Ali
  • 56,466
  • 29
  • 168
  • 265
  • Thanks a lot! This is very useful. The 'multistart' option is neat as I was doubting my initial values. I was very convinced about the model, so I tried a system of two equations. This sometimes gives less rubbish solutions like (w=1, x=4), but it is sensitive to changes in coefficients and the solution to w is always 1 for some reason. I'll get on to AMPL! Thanks again. – Divergent-Economist Mar 19 '14 at 14:26
  • @Divergent-Economist Your welcome; don't forget to chose an answer and accept it. And good luck! – Ali Mar 19 '14 at 15:05
0

This method will often rapidly find approximate solutions, NMinimize the sum of squares with constraints.

In[2]:= NMinimize[{
   (7*w - (7*w + 5*x + 2*y + z)*(0.76 + 0.12*Log[w] - 0.08*Log[x] - 
     0.03*Log[y] - 0.07*Log[7*w + 5*x + 2*y + z]))^2 +
   (5*x - (7*w + 5*x + 2*y + z)*(0.84 - 0.08*Log[w] + 0.11*Log[x] - 
     0.02*Log[y] - 0.08*Log[7*w + 5*x + 2*y + z]))^2 +
   (2*y - (7*w + 5*x + 2*y + z)*(-0.45 - 0.03*Log[w] - 0.02*Log[x] + 
     0.05*Log[y] + 0.12*Log[7*w + 5*x + 2*y + z]))^2 +
   (1*z - (7*w + 5*x + 2*y + z)*(-0.16 + 0*Log[w] + 
     0.03*Log[7*w + 5*x + 2*y + z]))^2, 
   w > 0 && x > 0 && y > 0 && z > 0}, {w, x, y, z}, 
   Method -> "RandomSearch"]

Out[2]= {9.34024*10^-12, {w->1.86998*10^-8, x->3.83383*10^-8, y->4.59973*10^-8, z->5.29581*10^-7}}
Bill
  • 3,664
  • 1
  • 12
  • 9
  • Thanks. Makes sense to look for global minimum using this! Though as Ali pointed out these solution values are not the ones I hoped to see. But at least your method is ensuring me that I haven't missed out on solutions to the system. – Divergent-Economist Mar 19 '14 at 14:35
  • I am always skeptical and I would have to think about that a lot longer before I could be convinced there were no other solutions than the one near zero. I don't know that the implementation details of RandomSearch are really documented so I still worry about this. But this technique seems to usually be fast and find a solution when other methods sometimes grind away until infinity. – Bill Mar 20 '14 at 04:41