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;
Improve this question
$\endgroup$

1 Answer 1

Reset to default
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*)
Improve this answer
$\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
    Commented 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
    Commented 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
    Commented 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
    Commented Jan 11, 2017 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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