0
$\begingroup$

I want to solve multiple equations so I thought using NMinimize. The basic idea is that if you have multiple equations with multiple variables :

$$\left\{\begin{split}&f(x,y,z)=0\\ &g(x,y,z)=0 \\& h(x,y,z)=0\\\end{split}\right.$$

You can find the solution by minimizing : $|f(x,y,z)|+|g(x,y,z)|+|h(x,y,z)|$.

So this is what I'm trying to do here in this code with the functions p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z],Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z], Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z].

My problem is that when I'm changing a little bit the condition : z>... the result becomes different.

Here is the introduction to the NMinimize :

a0=0.0748294;
b0=0.629316;
h[x_, y_] = 
  a*x*Log[x] + b*(1 - x - y)*Log[1 - x - y] + c*y*Log[y] - x^2 - d*x*y;
"------------------------------"
Hx[x_, y_] = D[h[x, y], x];
Hy[x_, y_] = D[h[x, y], x];
hxx[x_, y_] = D[h[x, y], x, x];
hxy[x_, y_] = D[h[x, y], x, y];
hyy[x_, y_] = D[h[x, y], y, y];
det[x_, y_] = Det[{{hxx[x, y], hxy[x, y]}, {hxy[x, y], hyy[x, y]}}];
p[x_, y_] = h[x, y] - x*Hx[x, y] - y*Hy[x, y];
h[a_, b_, c_, d_, x_, y_] = h[x, y];
p[a_, b_, c_, d_, x_, y_] = p[x, y];
Hx[a_, b_, c_, d_, x_, y_] = Hx[x, y];
Hy[a_, b_, c_, d_, x_, y_] = Hy[x, y];
det[a_, b_, c_, d_, x_, y_] = det[x, y];
x = 0.001;
y = 0.001;
w = 0.4;

And here are the several NMinimize changing a little bit the condition. As you can see, the condition I'm changing is z>0.2,z>0.3,z>0.4. And as you can see all the results give $z>0.4$ so none of them contradicts any of the conditions.

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.2 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.3 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.4 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

The results :

-> {0., {c -> 0.584662, d -> 2.09663, z -> 0.486347}}

-> {1.94289*10^-16, {c -> 0.560407, d -> 2.26232, z -> 0.507207}}

-> {1.8735*10^-16, {c -> 0.114809, d -> 3.31333, z -> 0.571615}}

As you can see they are all different. How can I do to find the exact result ?

EDIT :

What I ended up doing is to write :

sol = Table[
  NMinimize[{((p[a0, b0, c, d, x, y] - 
          p[a0, b0, c, d, w, z])^2 + (Hx[a0, b0, c, d, x, y] - 
          Hx[a0, b0, c, d, w, z])^2 + (Hy[a0, b0, c, d, x, y] - 
          Hy[a0, b0, c, d, w, z])^2)*1000000, 
    c > 0 && d > 0 && z > nn && z < 1 - w && 
     det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
   AccuracyGoal -> 20, PrecisionGoal -> 22, 
   MaxIterations -> 1000], {nn, 0.3, 0.575, 0.025}]

And I plotted the results, as one can see it converges :

ListPlot[sol[[All, 1]]]
ListPlot[Table[z /. sol[[n, 2, 3]], {n, 1, Length[sol]}]]
ListPlot[Table[d /. sol[[n, 2, 2]], {n, 1, Length[sol]}]]
ListPlot[Table[c /. sol[[n, 2, 1]], {n, 1, Length[sol]}]]

enter image description here

I'm not sure it's a good method, even now convergence is generally a good sign.

$\endgroup$
2
  • 1
    $\begingroup$ Usually, for differentiability purposes, one considers minimizing $(f^2+g^2+h^2)/2$, after which one can use something like Levenberg-Marquardt. Also (and apologies for the self-promotion), have you already seen this? $\endgroup$
    – J. M.'s torpor
    Apr 23 '20 at 13:12
  • $\begingroup$ @J.M. Thank you ! I just used your code FindAllCrossings3D but I get multiple results when I'd like to have only one and the value of Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]] is a bit higher than with NMinimize. So it's a great code but I don't manage to get the right results for my case. How can I do this Levenberg-Marquardt method with Mathematica ? $\endgroup$
    – J.A
    Apr 23 '20 at 14:06
2
$\begingroup$

Edit Corrected error, not taking second equation into account.

In this case solution with NMinimize is not good, since you get a whole curve of solutions, not only single points.

If i got it right the difference of the p and the Hx have to be zero. Hy is the same as Hx. I rationalized to get a solution with Solve.

sol = Solve[{p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z] == 0, 
Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z] == 0} // 
Rationalize[#, 0] &, {c, d}]

(csol[z_] = c /. First@sol) // N

(*   (7.99877*10^-29 (-1.16357*10^37 - 3.27364*10^37 z - 
 1.91963*10^37 Log[0.6\[VeryThinSpace]- 1. z] + 
 1.92156*10^37 z Log[0.6\[VeryThinSpace]- 1. z]))/(1.68711*10^7 + 
 2.44235*10^9 z Log[z])   *)

(dsol[z_] = d /. First@sol) // N

(*   -((9.51773*10^-11 (3.68704*10^12 + 
6.61204*10^12 Log[0.6\[VeryThinSpace]- 1. z]))/(-1. + 1000. z))   *)

.

ParametricPlot3D[{csol[z], dsol[z], z}, {z, 1/5, 3/5}, 
 AspectRatio -> 1, AxesLabel -> {c, d, z}, PlotStyle ->  Thick, 
 RegionFunction -> 
 Function[{c, d, 
z}, -d^2 - 
  157968671/(138459256 (3/5 - z)) + (157329 d)/(125000 (3/5 - 
       z)) - (3625853 c)/(2000000 z) + (157329 c)/(250000 (3/5 - 
       z) z) > 0 && c > 0]]

enter image description here

Test with the point minima you got

dsol[0.507207]
(*   2.26232   *)

csol[0.507207]
(*   0.560406   *)
$\endgroup$
1
  • $\begingroup$ I think I made a mistake in the first post in the definition but I get your method thank you ! $\endgroup$
    – J.A
    Apr 24 '20 at 9:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.