1
$\begingroup$

I'm trying to find the minimun value of Etot with respect to the variables bvalue and cvalue. bvalue has a range set by two constants bmin and bmax. But cmin, cmax are values with respect to bvalue.

How can I change this code to get the minimum value of Etot and the corresponding avalue, bvalue?

(* range of b *)
    Clear[a, a00, b, b00, c, c00, z0, \[Beta], v0, v1, Etot, Ea, Eb, Eg, \
    \[Gamma]lv, k, g, \[Rho], avalue, bvalue, bmin, bmax, cvalue, cmin, \
    cmax, z0value, z1, sol1, sol2, S1, S2, S3, r]
    Off[Reduce::ratnz];
    Off[Solve::ratnz];
    \[Beta] = 135000;
    v0 = 5*10^-9;
    \[Gamma]lv = 72.5*10^-3;
    k = 0.94;
    g = 9.8;
    \[Rho] = 1000;
    (*
    bvalue=1.09\[Times]10^-3;
     cvalue=0.82\[Times]10^-3; 
    *)
    sol1 = Reduce[{a/b^2 + 1/a == (2 b)/a^2 + 
           2/3 \[Beta] b && (3 v0)/(4 \[Pi] a^2) < 
          b < (3 v0)/(2 \[Pi] a^2)}, a, Reals];
    {a, {bmin, bmax}} = 
      N[sol1] /. {HoldPattern[
          And[_[b1_, ___, b2_], a == a0_]] :> {a0, {b1, b2}}};
    afn[b0_] := a /. b -> b0;
    avalue = afn[bvalue];
    v1fn[a_, b_] := 2/3 \[Pi] a^2 b;
    v1 = v1fn[avalue, bvalue];

(* range of c *)
sol2 = Part[
   Solve[v0 == 
     2/3 \[Pi] avalue^2 bvalue + 2/3 \[Pi] avalue^2 cvalue - 
      1/3 \[Pi] avalue^2/cvalue^2 (2 cvalue^3 - 3 cvalue^2 z0 + z0^3),
     z0], 2];
{cmin, cmax} = {3/(2 \[Pi] avalue^2) (v0 - v1), bvalue};
z1 = z0 /. sol2;
z0fn[c0_] := z1 /. c -> c0
z0value = z0fn[cvalue];

(*Gibbs evergy from here *)
r[a_, c_, z0_] := a (1 - z0^2/c^2)^(1/2);
e[x_, y_] := (1 - (y/x)^2)^(1/2)
\[Eta][x_, y_, z_] := ((e[x, y] x z)^2 + y^4)^(1/2)

S[x_, y_, z_] := \[Pi]/(
  y^2 e[x, y]) (y^2 (e[x, y] x^2 + y^2 Log[e[x, y] x y + x y]) - 
    e[x, y] x z \[Eta][x, y, z] - 
    y^4 Log[e[x, y] x z + \[Eta][x, y, z]])

S1 = S[avalue, bvalue, 0];
S2 = S[avalue, cvalue, 0];
S3 = S[avalue, cvalue, z0value];

(* the energies and the values of a,b,c,z0 *)
Ea = \[Gamma]lv (S1 + S2 - S3)
Eb = \[Gamma]lv \[Pi] r[avalue, cvalue, z0value]^2 k
Eg = \[Rho] g (2/3 \[Pi] avalue^2 bvalue (3/8 bvalue + z0value) + 
    1/12 \[Pi] avalue^2/bvalue^2 z0value^2 (6 cvalue^2 - z0value^2))
Etot = Ea + Eb + Eg

{avalue, bvalue, cvalue, z0value};
{bmin, bmax, cmin, cmax};

(* Energy Minimization*)

f[b_?NumericQ,c_?NumericQ]:=Etot/.{bvalue\[Rule]b,cvalue\[Rule]c};
cm[b_]:=cmin/.bvalue\[Rule]b;
cmx[b_]:=cmax/.bvalue\[Rule]b;
f[0.00109,0.00082]

NMinimize[{Re[f[bvalue,cvalue]],cm[bvalue]<cvalue<cmx[bvalue]},{{\
bvalue,bmin,bmax},cvalue}]

The problem is that the output of f[0.00109,0.00082] gives me a complex number when it should be same with when I put in the commented out avalue, bvalue(towards the top of the code when I define variables.

The output of Minimize gives me {-0.0597908, {bvalue -> 0.00122385, cvalue -> 2.16242*10^-11}} . The expected answer is something close to;

bvalue=1.09\[Times]10^-3;
cvalue=0.82\[Times]10^-3;
Etot=1.06856*10^-6;
$\endgroup$

1 Answer 1

2
$\begingroup$

NMinimize seems to be having a hard time dealing with the constrain equation on cvalue. One possible solution modifies the objective function such that minimizing necessarily satisfies the required constraints. Defining such a function

consteq[bvalue_, cvalue_] := If[cm[bvalue] < cvalue < cmx[bvalue], 0, 10]

Now NMinimize takes the form

NMinimize[{Re[f[bvalue, cvalue]] + consteq[bvalue, cvalue], 
  bmin < bvalue < bmax}, {bvalue, cvalue}]
(*{-1.92845*10^-6, {bvalue -> 0.00127286, cvalue -> 0.0000942478}}*)

Changing the Method in NMinimize to DifferentialEvolution improves the minimization solution a bit

NMinimize[{Re[f[bvalue, cvalue]] + consteq[bvalue, cvalue], 
  bmin < bvalue < bmax}, {bvalue, cvalue}, 
 Method -> "DifferentialEvolution"]
(*{-1.96088*10^-6, {bvalue -> 0.00128478, cvalue -> 0.0000551245}}*)

You mention the following at the end of your post:

bvalue=1.09*10^-3;
cvalue=0.82*10^-3;
Etot=1.06856*10^-6;

However, checking with the definitions provided in the question does not yield the expected result

Etot /. {bvalue -> 1.09*10^-3, cvalue -> 0.82*10^-3}
(*-1.4596*10^-6 + 0. I*)
$\endgroup$
4
  • $\begingroup$ Thank you for the reply. I tried Etot /. {bvalue -> 1.09*10^-3, cvalue -> 0.82*10^-3} too and got the same result as you too(-1.4596*10^-6 + 0. I). I don't know why i'm getting this complex number. Because if i define the values of avalue, bvalue at the top of the code(the part that I have commented out) Etot returns 1.06856*10^-6 just fine $\endgroup$
    – Jun
    Jan 11, 2017 at 0:16
  • $\begingroup$ Also I tried inputting bvalue and cvalue by your method /. {bvalue -> 1.09*10^-3, cvalue -> 0.82*10^-3} and found out that S3 and r give different outputs from when you define them at the beginning. S1 and S2 were the same. So i suspect z0 to be the source of the problem but can't figure it out.. $\endgroup$
    – Jun
    Jan 11, 2017 at 0:33
  • $\begingroup$ At the first comment It should be; Because if i define the values of 'bvalue,cvalue' at the top of the code(the part that I have commented out) Etot returns 1.06856*10^-6 just fine. Can't edit it for some reason $\endgroup$
    – Jun
    Jan 11, 2017 at 0:34
  • $\begingroup$ omg I finally figured it out! I took the second element of sol2 for z0 because that was the right value for when I define bvalue,cvalue at the beginning of the code but I tried all three roots and found that the third one was the right one. Changing sol2=Part[..,3] gives me the answer I wanted awesome. $\endgroup$
    – Jun
    Jan 11, 2017 at 1:03

This site is temporarily in read-only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .