0
$\begingroup$

I'm trying to write a function of a variable as the argument of a minimization where that variable comes to play.

Here is the code :

f[x_, a_, b_] = x*Log[x] + (1 - x)*Log[1 - x ] - a*x^2 + b*x^4;
mu[x_, a_, b_] = D[f[x, a, b], x];
h[x_, a_, b_] = D[f[x, a, b], x, x];
solh[a_, b_] := {Min[x /. NSolve[{h[x, a, b] == 0, 0 < x < 1}, x]], 
   Max[x /. NSolve[{h[x, a, b] == 0, 0 < x < 1}, x]]};

p[x_, a_, b_] = f[x, a, b] - x*mu[x, a, b];

fm[a_, b_] := 
  NMinimize[{(mu[x0, a, b] - mu[x1, a, b])^2 + (p[x0, a, b] - 
        p[x1, a, b])^2, 
    0.0001 < x0 < solh[a, b][[1]] && 
     solh[a, b][[2]] < x1 < 0.9999}, {x0, x1}];

amin = 2;
amax = 82;
da = 5;
bmin = 10;
bmax = 100;
db = 20;

X = Table[Table[fm[a, b], {a, amin, amax, da}], {b, bmin, bmax, db}];
A = Table[a, {a, amin, amax, da}];
X0 = Transpose[{A, x0 /. X[[All, 2]]}];
X1 = Transpose[{A, x1 /. X[[All, 2]]}];
ListPlot[{X0, X1}]

But I'm then getting the error :

NMinimize::bcons: The following constraints are not valid
$\{0.0001<x0,x<x1,x0<x,x1<0.9999\}. $ Constraints should be equalities, inequalities, or domain specifications involving the variables.

So I was wondering how to proceed !

Thanks in advance

$\endgroup$
1
$\begingroup$

For some values of a and b there is no solutions in the range 0<x<1, so solh[a_, b_] shows {x, x}, hence further errors. You can change this to:

solh[a_, b_] := 
  With[{sol = NSolve[{h[x, a, b] == 0, 0 < x < 1}, x][[;; , 1, 2]]}, 
   If[sol != {}, MinMax[sol], {0.00011, 0.99989}]];

Now in case of no solutions due to restrictions on x the output will be {0.00011, 0.99989}, you can tune these values. Then we can simplify your Tables:

points = Flatten[
  Table[{a, x0, x1} /. fm[a, b][[2]], {b, bmin, bmax, db}, {a, amin, 
    amax, da}], 1];

and show ListPlot:

ListPlot[{points[[;; , {1, 2}]], points[[;; , {1, 3}]]}, PlotLegends -> {x0, x1}]

enter image description here

$\endgroup$

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.