1
$\begingroup$

I am trying to find a set of parameters $\{\epsilon,\epsilon_b\}$, so that the two local minima have exactly the difference of 15.1 (see code below) and the left sided local minimum and the local maximum have a distance of 2.5 (see code below). I am able to solve this problem in form of two ContourPlots and Show, but I am not able to find the intercept of the two ContourLines with NSolve. Maybe someone can help me with that and explain why I can solve this problem graphically, but not numerically.

The code:

U = Function[{r, ϵ, ϵb},E^(-((2(r-2σ)^2)/[Sigma]^2)) ϵb + 4 ϵ (-(σ^6/r^6) + σ^12/r^12) + (lc Log[(lc (1 - r^2/lc^2) Sinh[(r (3 - r^2/lc^2))/(lc (1 - r^2/lc^2))])/(r (3 - r^2/lc^2))])/(β σ)];

T = 298.15;
kB = 8.3145*10^-3;
β = 1/(kB T);
lc = 65;
σ = 2.5;

xfCoord[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}]
xf[ϵ_?NumberQ, ϵb_?NumberQ] := r /. FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}][[2]]
Uf[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}][[1]]
xuCoord[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMinimum[{U[r, ϵ, ϵb],3 σ < r < lc}, {r,3 σ}]
xu[ϵ_?NumberQ, ϵb_?NumberQ] :=r /. FindMinimum[{U[r, ϵ, ϵb],3 σ < r < lc}, {r, 3 σ}][[2]]
Uu[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb], 3 σ < r < lc}, {r,3 σ}][[1]]
xbCoord[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMaximum[{U[r,ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}]
xb[ϵ_?NumberQ, ϵb_?NumberQ] :=r /. FindMaximum[{U[r,ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}][[2]]
Ub[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMaximum[{U[r, ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}][[1]]

Δx[ϵ_?NumberQ, ϵb_?NumberQ] := xb[ϵ, ϵb] - xf[ϵ, ϵb]
ΔG[ϵ_?NumberQ, ϵb_?NumberQ] := Uu[ϵ, ϵb] - Uf[ϵ, ϵb]

p1 = ContourPlot[Δx[ϵ, ϵb] == 2.5, {ϵ, 2, 20}, {ϵb, 1, 20},FrameLabel -> {"ϵ","ϵb"}];
p2 = ContourPlot[ΔG[ϵ, ϵb] == 15.1, {ϵ, 2, 20}, {ϵb, 1, 20},FrameLabel -> {"ϵ","ϵb"},ContourStyle -> Green];
Show[p1, p2]

ContourLines

This now does not work:

NSolve[{Δx[ϵ, ϵb] == 2.5, ΔG[ϵ, ϵb] == 15.1}, {ϵ, ϵb}]

It just gives back the same expression, meaning it can not evaluate it, I guess.

$\endgroup$
3
  • $\begingroup$ I'm not sure you gave us all your code. I'm seeing [Sigma] in your definition for U. What is [Sigma]? $\endgroup$ Oct 14, 2015 at 19:06
  • $\begingroup$ $\sigma = 2.5$ (defined in the second block of the code above). $\endgroup$ Oct 14, 2015 at 19:53
  • $\begingroup$ Then the slash in front of [sigma] should be \ instead of /. So \[sigma]. $\endgroup$ Oct 14, 2015 at 20:14

2 Answers 2

2
$\begingroup$

That's because NSolve[ ] isn't the right tool for the job as it deals primarily with linear and polynomial equations.

Try FindRoot[ ] instead:

FindRoot[{Δx[ϵ, ϵb] == 2.5, ΔG[ϵ, ϵb] == 15.1}, {{ϵ, 5}, {ϵb, 5}}]

(* {ϵ -> 14.6863, ϵb -> 5.28574} *)
$\endgroup$
1
$\begingroup$
U = Function[{r, ϵ, ϵb}, 
   E^(-((2 (r - 2 σ)^2)/σ^2)) ϵb + 
    4 ϵ (-(σ^6/r^6) + σ^12/
        r^12) + (lc Log[(lc (1 - 
             r^2/lc^2) Sinh[(r (3 - r^2/lc^2))/(lc (1 - r^2/lc^2))])/(r (3 - 
             r^2/lc^2))])/(β σ)];

T = 298.15;
kB = 8.3145*10^-3;
β = 1/(kB T);
lc = 65;
σ = 2.5;

To restrict a function's arguments to numeric values use NumericQ rather than NumberQ. NumberQ will not accept numeric symbols (e.g., π, E, GoldenRatio).

NumberQ /@ {π, E, GoldenRatio}

(*  {False, False, False}  *)

xfCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMinimum[
   {U[r, ϵ, ϵb], 0 < r < 2 σ},
   {r, σ}];
xf[ϵ_?NumericQ, ϵb_?NumericQ] :=
  r /. FindMinimum[{U[r, ϵ, ϵb], 
      0 < r < 2 σ}, {r, σ}][[2]];
Uf[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMinimum[
    {U[r, ϵ, ϵb], 0 < r < 2 σ},
    {r, σ}][[1]];
xuCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMinimum[{U[r, ϵ, ϵb], 3 σ < r < lc},
   {r, 3 σ}];
xu[ϵ_?NumericQ, ϵb_?NumericQ] :=
  r /. FindMinimum[
     {U[r, ϵ, ϵb], 3 σ < r < lc},
     {r, 3 σ}][[2]];
Uu[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMinimum[
    {U[r, ϵ, ϵb], 3 σ < r < lc},
    {r, 3 σ}][[1]];
xbCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMaximum[
   {U[r, ϵ, ϵb], σ < r < 0.9 lc},
   {r, 2 σ}];
xb[ϵ_?NumericQ, ϵb_?NumericQ] :=
  r /. FindMaximum[
     {U[r, ϵ, ϵb], σ < r < 0.9 lc},
     {r, 2 σ}][[2]];
Ub[ϵ_?NumericQ, ϵb_?NumericQ] :=
  FindMaximum[
    {U[r, ϵ, ϵb], σ < r < 0.9 lc},
    {r, 2 σ}][[1]];

Δx[ϵ_?NumericQ, ϵb_?NumericQ] :=
  xb[ϵ, ϵb] - xf[ϵ, ϵb];
ΔG[ϵ_?NumericQ, ϵb_?NumericQ] :=
  Uu[ϵ, ϵb] - Uf[ϵ, ϵb];

ContourPlot[
 {Δx[ϵ, ϵb] == 
   2.5, ΔG[ϵ, ϵb] == 15.1},
 {ϵ, 2, 20}, {ϵb, 1, 20},
 FrameLabel -> {"ϵ", "ϵb"}]

enter image description here

Use FindRoot. The initial estimates can be taken from the plot

FindRoot[{Δx[ϵ, ϵb] == 
   2.5, ΔG[ϵ, ϵb] == 15.1},
 {{ϵ, 15}, {ϵb, 5}}]

(*  {ϵ -> 14.6863, ϵb -> 5.28574}  *)
$\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.