(* 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]
\[Beta] = 135000;
v0 = 5*10^-9;
\[Gamma]lv = 72.5*10^-3;
k = 0.94;
g = 9.8;
\[Rho] = 1000;
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*)



So 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?

(I have changed the brackets in S function as you mentioned and the variable r into a function. Because when I tried putting in the commented out values of bvalue and cvalue(towards the top part of the code where I define my variables), Etot was a complex number. I tried changing it into a function and its location and now Eb and Etot are both positive numbers as should be.

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 the avalue, bvalue.

Also the output of Minimize gives me {-0.0597908, {bvalue -> 0.00122385, cvalue -> 2.16242*10^-11}} When it should give me sth close to

  • $\begingroup$ How do you change them like that? The code I wrote on my mathematica were like those too but when I copy and pasted them here they changed automatically like the ones before. $\endgroup$
    – Jun
    Jan 9 '17 at 15:22
  • 1
    $\begingroup$ Chrome Extension Have you seen this? $\endgroup$
    – Feyre
    Jan 9 '17 at 15:24
  • $\begingroup$ I only read the second answer. The first one was too hard for me with all the stuff like Options[#, {Cubics, Quartics}] & /@ {Reduce, Solve} I don't know what &/@ means $\endgroup$
    – Jun
    Jan 9 '17 at 15:33
  • $\begingroup$ I downloaded the extension but the buttons don't show.. $\endgroup$
    – Jun
    Jan 10 '17 at 8:38

change the last few lines to look like this:

Etot = Ea + Eb + Eg;
f[b_?NumericQ, c_?NumericQ] := 
  Etot[[1, 1]]  /. {bvalue -> b, cvalue -> c};
cm[b_] := cmin /. bvalue -> b;
cmx[b_] := cmax /. bvalue -> b;
NMinimize[{Re[f[bvalue, cvalue]], 
  cm[bvalue] < cvalue < cmx[bvalue]}, {{bvalue, bmin, bmax}, cvalue}]

notice Etot is a 2d list so I just took the first part. You probably want to review your code to determine why its a list in the first place. (Do you really want the curly brackets in the S function?) It is also complex with a small imaginary part so I just took the real part.

  • $\begingroup$ Thank you for replying. But i'm still having trouble. I edited my question. The commented out values of avalue, bvalue and the output of Etot from when you put those in are what i'm expecting as the answer. Etot is suppoed to be a positive real number. $\endgroup$
    – Jun
    Jan 10 '17 at 8:18
  • $\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 '17 at 1:04

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.