在 Mathematica 中求解四个方程组（带对数）

Question

我正在尝试求解一个包含四个变量的四个方程组。我已经阅读了许多关于类似问题的主题，并尝试遵循这些建议。但我认为这有点乱，因为这里的日志和交叉产品。这是确切的系统：

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

（仅供参考-我是一名年轻的经济学家，这是消费者需求系统的延伸。）理论上，我们知道这个系统存在一个独特的解决方案，它是积极的。

试试

• Solve & NSolve ：因为应该有解决方案，所以我尝试了这些，但都不起作用。我猜系统有太多日志需要处理。

• FindRoot ：我从数据中获得的初始值 (14,15,10,100) 开始。 FindRoot 返回最后一个值（不满足我的系统）和以下消息。

  FindRoot::lstol：线搜索将步长减小到 AccuracyGoal 和 PrecisionGoal 指定的容差范围内，但无法.....

我尝试了不同的初始值，包括 FindRoot 返回的值。我试图在每一步分析解决方案值的模式。我没有看到任何模式，但注意到 z 值在此过程的早期变为负数。所以我对这些值设置了界限。这只是在最小值 0.1 处停止代码。我还尝试了指数系统而不是日志，同样的问题。

 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 ：因为我们可以将此问题写为最小化问题，所以我尝试了这个（按照这里的建议）。返回的值没有将系统或方程收敛到零。我只尝试了前两个方程，并且那种收敛到零。

希望这对这里的专家来说足够有吸引力！有什么想法我应该如何找到解决方案，或者为什么我不能？这是我第一次使用 Mathematica，不幸的是我第一次凭经验解决系统/优化问题！非常感谢。

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

Answer 1

我无法访问 Mathematica，我将您的方程放入了AMPL，它对学生免费。这是我所做的：

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

如果我只要求变量是非负的，例如，我会得到垃圾

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

或

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

从数值上看，这些确实是解决方案，它们以相当高的精度满足方程，尽管我很确定这不是您要寻找的。这表明您的模型存在问题。

如果我稍微增加变量的下限：

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

我得到，即使是多启动，Ipopt 3.11.6: Converged to a locally infeasible point. Problem may be infeasible. 这再次表明您的模型方程存在问题。

恐怕你将不得不修改你的模型。

---

这不会解决您的模型方程的问题，但我会引入新变量：a=log(w), b=log(x), c=log(y), d=log(z)。然后w=exp(a) 等等。它的优点是函数评估不会由于对数函数的负参数而失败。

我可能还会为 (7*w+5*x+2*y+z) 引入一个新变量，只是为了使方程更紧凑。

这些新变量都不会解决您的模型方程的上述问题。

---

如果您真的是第一次使用 Mathematica，那么最好使用 AMPL 和 IPOPT；这些工具是为求解方程和优化问题而定制的。如果您有疑问而不是 Stackoverflow，我建议您使用 AMPL mailing list；您将在邮件列表中获得更好的答案。

Answer 2

这种方法通常会快速找到近似解，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}}