# NMinimize : different results when changing a little bit the constraints

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


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

• 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? Apr 23 '20 at 13:12
• @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 ?
– J.A
Apr 23 '20 at 14:06

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


Test with the point minima you got

dsol[0.507207]
(*   2.26232   *)

csol[0.507207]
(*   0.560406   *)

• I think I made a mistake in the first post in the definition but I get your method thank you !
– J.A
Apr 24 '20 at 9:21